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ABSTRACT 



System identification concerns the mathematical. modeling of a system based upon 
its input and output. It allows the development of a mathematical description when all 
that is available is the result of a process or the output of a system and not the process 
or system itself. 

The purpose of this thesis is to develop algorithms for modeling systems as 
autoregressive-moving-average processes using the method of instrumental variables, a 
modification of the ordinary least-squares technique, and a multichannel method based 
upon processing the input and output data by separate infmite-impulse-response filters. 
The methods developed are tested by computer simulation using several second and 
third-order test cases and the results are presented. 
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1. INTRODUCTION 



A. SYSTEM IDENTIFICATION BASICS 

System identification concerns the modeling of systems as sets of mathematical 
equations based upon the input and output of the system. [Ref 1: pp. 3-6]. It allows 
a model to be developed when all that is available is the result of a process or the output 
of a system and not the process or the system itself System identification is an impor- 
tant area of study. Solution of the modeling problem offers many alternatives for the 
continued study of the system. Among these are: 

• Nondestructive analysis of the system. 

• Simulation studies using the model. 

• Easy adaptation of the model to changing system environment. 

• Spectral analysis of the system. 

Modeling can simulate the system's operation at a fraction of the cost of actual 
system operation. Complex operations not possible with the actual system for fear of 
damaging it or personal injury can be simulated. This can expose how the system will 
operate in adverse conditions not normally experienced. In speech processing, modeling 
the speech process has the potential for significantly reducing the amount of information 
necessary to store in order to reproduce the speech. 

The modeling process shown in Figure 1 on page 2 assumes the unknown system's 
input and the output data are available for processing. In many cases, if the system's 
input is unknown or data is not available, a white noise input can be used in its place. 
The modeling process uses the input and output data to find a set of parameters which 
closely approximate the operation of the system. The better the identification technique, 
the more closely the model follows the performance of the actual system. 

Many types of models are available. This thesis investigates a linear parametric 
model that can be described by difference equations. This type of model lends itself well 
to simulation on a digital computer. The frequency characteristics of the system deter- 
mined from the parameters of these types of models are more accurate than what can 
be determined from classical means such as FFTs. This is because classical methods use 
windows which assume data beyond their extent is zero [Ref. 2: p. 173). This is not a 
realistic assumption. Models in this category include the moving-average (MA) model, 
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Figure 1. System identification problem 

the autoregressive (AR) model, and the autoregressive-moving-average (ARMA) model. 
In the frequency domain, MA processes are characterized by sharp nulls and smooth 
peaks and AR processes are characterized by smooth nulls and sharp peaks. ARMA 
processes have sharp peaks and sharp nulls [Ref. 2: p. 173]. An advantage of the MA 
process is its inherent stability. .An advantage of the .AR process is the large number of 
algorithms already available for modeling systems. An advantage of the ARMA process 
is that it uses far fewer parameters than either the M.A or AR process alone to model a 
system. This satisfies the general requirement to reduce the complexity of the model. 

In addition to a large variety of models, there are two processing modes: block and 
sequential. 

Block processing uses a fixed length block of data in the parameter estimation 
process. It ignores data before and after the block. This is not a real time processing 
method because all data must be available before processing can start. Block processing 
generally involves inversions of data matrices whose sizes are on the order of 
(.V-t- M) X (.V + }[) where .V is the order of the .AR process and M is the order of the 
MA process. 



Sequential processing uses new data to update the parameter estimations. It starts 
by initializing an estimate of the inverse of the data covariance as as a diagonal matrix. 
It uses each new data point to update this matrix. Then it updates the parameter esti- 
mates using the updated inverse data covariance matrix. It is a real time method. I'he 
algorithm to implement the sequential processing method is generally more complex 
than the block method but less computationally intensive because the matrix inversions 
are not required. 

This thesis concerns only systems represented by discrete time data uniformly sam- 
pled at a sufficient rate to meet the Xyquist criteria. 

The work in this thesis assumes that the input data is a wide-sense stationarx' ran- 
dom sequence. Tests of the algorithms used a pseudorandom Gaussian input with a 
mean of zero and a variance of one. 

B. PROBLEM STATEMENT 

The purpose of this thesis is to develop algorithms for modeling systems as ARMA 
processes using the method of instrumental variables (IV) and a multichannel approach. 
Tests of the methods will be conducted to determine the accuracy of their results and the 
speed with which they converge. 

The IV approach is a modification of the method of ordinary' least squares. This 
approach is developed first as a block processing case and then converted to a sequential 
processing case. Tests are conducted of only the sequential processing case. 

Using a multichannel scheme allows the input and output data of the unknown 
system to be processed separately. This reduces the sizes of the data matrices involved 
in the modeling process. Both block and sequential processing cases are formulated but 
only the block processing case is tested. 

C. OVERVIEW OF THESIS 

Chapter 2 is about AR.MA modeling. It also presents a detailed derivation of the 
method of ordinarx’ least squares because it forms the basis on which other modeling 
techniques depend. 

Chapter 3 presents a modified least-squares approach called the method of instru- 
mental variables. It is attractive due to its simplicity and good noise performance. 
Chapter 3 presents results of using this method on sex’eral second and third-order test 
systems. 
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Chapter 4 presents a new multichannel approach to ARMA modeling. This ap- 
proach is presented in block and sequential processing forms. This chapter also presents 
several adaptations of the block form which improve its speed of convergence. 

Chapter 5 contains a summaiy of the thesis and lists topics for further research. 

The appendix contains the programs used to test the sequential IV algorithm and 
the block multichannel iterative algorithm. Subroutines common to both programs are 
grouped together and listed at the end of the appendix. 
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II. ARMA MODELING 



A. ARMA PROCESSES 

Modeling as an autoregressive-moving-average (ARMA) process has the potential 
for achieving a close fit to the system using a reduced order over that which a moving 
average or an autoregressive model alone could achieve. ARMA modeling is concerned 
with finding a set of AR parameters and .MA parameters which combined describe an 
ARMA process that approximates the characteristics of a target system. 

The general form of the ARMA model is shown in Figure 2 on page 6. The output 
at time n, y(n), is a linear combination of past outputs and past and present inputs. 
The a, and b, are constants referred to as tap weights. The parameters form the MA 
part of the ARM.A. model. The b, parameters form the AR part. In equation form the 
output of the ARMA system is represented by the following difference equation; 



;V M 

y{n) = - Y^biy{n - t) + u{ri - /) ( 2 . 1 ) 

i=l (=0 

where A' is the order of the AR part of the ARMA model and M is the order of the MA 
part of the ARM.A model. This means the ARMA output at the current time depends 
on the last .V values of the ARMA output. The N weighting parameters determine 
exactly how the new output depends on the past outputs. The M a, weighting parame- 
ters determine how the new output depends on the current and M — 1 past inputs. 

Equation (2.1) in vector form becomes; 

y = x^e ( 2 . 2 ) 

where x is a (A' -I- A/ 4- 1) x 1 vector of input and output data values given by; 

x = i -y{n - I) -y{n - 2) ... -y(a - A") x(«) x(« - 1) ... x{n - A/)]^ (2.3) 

and 0 is a (A’ -I- M -I- 1) x 1 vector of the AR and MA tap weights given by; 

9 = Lb, b, ... b^^ - %/]"' ( 2 - 4 ) 
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For iV + L — I data points available for y and M + L data points available for u, we 
can write a block equation which gives the value of the output at progressive sampling 
times: 



6 



><« - L + 1) 

y{n - L + 2) 



>■(«) 



The row in equation (2.5) is the value ofj- at time « + / based on output data available 
through time n + i — 1 and input data available through n + /. The row is identically 
equation (2.2) at time n + i. In vector form equation (2.5) becomes; 

y = XO (2.6) 

where 6 is defined in equation (2.4): y, the vector of output values, is given by; 

y = lv{n -L+l) y{n - L + 2} ... y(«)]^ (2.7) 

and X is a partitioned matrix with rows comprised of data vectors exactly like equation 
(2.3) only shifted in time. At successive sampling times, when new data is obtained, data 
used to calculate the previous output shifts one column to the right. The new data fills 
in the left mosty and u columns. The matrix X is given by; 





-y{n-L) ... 


+1) 


u{n — A + 1) 


u{n — /I +1) 




— T + 1) 


-y{n - rj +2) 


u{n — L + 2) 


u{n — n +2) 


x = 


-T(« - 1) - 


-y{n - A') 


u{n) 


u{n — M) 


where t] is defined as + L and 


M is defined as M + L. 







\\n+ 1), 
x\n + 2) 




• 




x\n + L) 



bx 






(2.5) 
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If the and b, are estimates of the true values of the AR and MA parameters, then 
the filter output will be an estimate of the true output. We use a hat over a variable (for 
example,/-) to indicate an estimated v'alue. Rewriting equation (2.6) using the estimated 
ARMA parameters results in; 



y = X0 (2.9) 

where Q is defined as: 

e = {b, b, ... bf, So Si ... S^]^ (2.10) 

and y is now the vector of estimated output values and is given by: 

y = [y(«-L+l) /(«-L + 2) ... y{n)f (2.11) 

Up until now, we have discussed estimating the output of a system given its input, 
past output, and an estimate of the parameters which describe it. If, however, we know 
the output and input of the system, based on these equations, we can use them to gen- 

A 

erate a set of S, and b, which produces an ARMA output which is the best possible es- 

A 

timate of the system output. Then the 4 and b, will be optimal parameters for describing 
the operation of the unknown system as an ARMA process. 

B. METHOD OF ORDINARY LEAST-SQUARES 

In this thesis we use the method of ordinary least- squares as the means of finding 
the optimum set of ARMA parameters. It is a well known modeling technique. It offers 
the advantage of being widely used in the scientific community for a variety of modeling 
problems. It has been applied successfully to a large number of modeling problems with 
good results and has been successfully applied to classes of problems for which other 
methods have failed. [Ref. 3: p. 4] 

To apply the method of ordinary least-squares to system identification we form the 
error between the actual system output and the estimated output generated by the 
ARMA model. This error is given by: 

e = y-y = y-Xd (2.12) 

where y is the vector of the actual system outputs given by equation (2.7) and y is the 
vector of ARMA outputs given by equation (2.11). [Ref 1: p. 176] 
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We then let the sum of the squares of the errors at the instances of time the meas- 
urements of the data were taken become a measure of how well the estimates approxi- 
mate the true system outputs. This measure of performance, or cost function, is denoted 
J . It is written in equation form as: 










(2.13) 



Replacing the error in equation (2.13) with its equivalent expression from equation 
(2.12) results in: 



7=y'y + O^XX'0-2O'Xy (2.14) 

Equation (2.14) shows that the performance measure is a function of the estimated pa- 
rameters. The criterion is to minimize the measure of performance by taking its deriva- 
tive with respect to the parameter estimates and setting it equal to zero. Then equation 
(2.14) becomes: 



= 0 = 0 + 2XX^0 - 2X^y (2.15) 

cO 

A 

Solving for 0, the parameters, gives us the result: 

0 = (X^X)-’x^y (2.16) 

Equation (2.16) is the ordinary least-squares solution for the optimum AR.MA pa- 
rameters. It provides the best possible description, in a least-squares sense, of the data 
source. The resulting parameters provide the closest fit to the actual input and output 
data of the system in the sense of least-squares errors. 

Equation (2.16) uses a block processing approach. The product of XTC must be 

A 

formed and then inverted in order to calculate 6 . In addition to being computationally 
intensive, the estimate cannot be updated when new data becomes available without re- 
calculating (X^)~* . A sequential update which does not require (X^)"‘ to be recalcu- 
lated is presented in the next chapter in the context of the instrumental variable method 
of least-squares. 
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III. INSTRUMENTAL VARIABLE METHOD OF SYSTEM 

IDENTIFICATION 



A. INTRODUCTION 

The instrumental variable (IV) method of system identification is a variation of the 
method of ordinary least-squares. Its attraction over ordinary least-squares is that there 
is no bias in estimating the parameters when dealing with noise [Ref 4: p. 406]. Also, 
this method is known to yield consistent estimates and remains as easy to use as the 
method of ordinary’ least-squares [Ref 3; p. 119). 

When an additive noise term is present in the observable output, y{n), the output is 
given by; 



>•(«) = w(/r) -i- v(rt) (3.1) 

Here w{n) represents the actual output of the system and r(«) represents the noise. When 
this noise has a non-zero mean, using the noise corrupted output to model the unknown 
system by the ordinary least-squares approach leads to inaccurate estimates of its pa- 
rameters. The parameters are referred to as biased estimates. [Ref 3: p. 119, Ref 1: pp. 
192-193, and Ref 5; p. 704). 

The IV method shown in Figure 3 on page 1 1 generates an estimate of the un- 
known system's output by processing the input data through an auxiliary model which 
closely approximates the unknown system. In our implementation of the IV method, 
the au.\iliary model is an AR.MA model. Its output is free of the noise affecting the 
unknown system. The IV method uses the auxiliarx’ model output (estimate), w, to cal- 
culate the parameters of the unknown system. Therefore the IV parameter estimates are 
not biased like those generated by the method of ordinaiy least-squares. 

The IV method assumes the existence of a matrix Z composed of the auxiliary^ 
model's input and output data which has the following two properties [Ref 4: p. 406]: 



lim ^ Z^s = 0 
,V->oo A' 


(3.2) 


lim -|rZ^X = Q 

W-^oo A 


(3.3) 



where s is the error in fitting the parameter estimates to the data and is given by: 
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v(n) 




Figure 3. Modeling by the instrumental variable method 



s = \-X6,y (3.4) 

and Q is a nonsingular square matrix. 

The first property means Z is orthogonal to the error. This leads to the cancellation 
of the bias term inherent in ordinary least-squares techniques [Ref. 4: p. 406]. The sec- 
ond property ensures the inverse of e.xists. Z is assumed to have the same structure 
and size as the data matrix X in equation (2.8). Its contents differ in that the noise 
corrupted system output >•(«) in X is replaced by the output of the auxiliary model fv(?z) 
in Z. The new data matrix Z is given by: 





1 

1 


— w{n — -fl) 


u(n-L+\) ... 


u{n — ^ -M) 




— w{n — L + 1) 


— w{n — t} -t-2) 


u{n — L + 2) 


u{n — -1-2) 


z = 


— \v{n — 1) 


— vv(/j — l\) 


u{n) 


u{n - .\f) 
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where rj is defined as .V+ L and n is defined as M + L. Comparing X in equation (2.8) 
and Z in equation (3.5), we note the substitution of w(a/) for y{n). Thus, we are now 
using estimates of the true output \v{n) instead of noise corrupted samples >•(«). 

To incorporate Z into the parameter estimation process we begin with equation 
(2.12), which we rewrite as: 



y = Xe + £ (3.6) 

A 

This equation says that the estimates of the output, given by X6, differ from the actual 
outputs, y, by some fitting error e. Multiplying equation (3.6) by Z^ yields: 

Z^y = Z^Xd + Ze (3.7) 

A 

Equation (3.3j ensures that Z^ can be inverted. Solving for 6 results in: 

fl = (Z^X)~’Z^y-(Z^X)“*Z^£ (3.8) 

The (Z^)"'Z^y term in equation (3.8) is the IV estimate of the parameters. It is written 
as: 



0,^,= (Z^X)-’Z^y (3.9) 

The (Z^)"'Z^£ term in equation (3.8) represents a potential bias in the estimate. The 
first property of the Z matrix, given in equation (3.2), ensures this bias goes to zero, 
asymptotically. Applying this property, equation (3.8) can be rewritten as: 

6 = (Z"'X)-'z"'y = 0;^. (3.10) 

Equation (3.10) gives an unbiased estimate of the ARM A parameters. [Ref 1: pp. 
192-193]: 

Other least-squares methods avoid the bias inherent in ordinary least-squares but 
they are more complicated than the IV method to implement [Ref 3: p. 119]. Although 
this thesis does not attempt an analysis of the IV method in the presence of noise, any 
practical system identification technique must deal with noise. Hence, the attraction of 
and the desire to use the IV method. 

Equation (3.10) represents the block processing case. It assumes A' -f L — 1 output 
samples and M + L input samples arc available. These samples are used to calculate an 
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estimate of the parameters. Samples beyond this range are not included in the esti- 
mation process. Block processing involves multiplication of two L x (A'-i- M + 1) ma- 
trices to form a third matrix. Then this third matrix must be inverted. This is a 
computationally intensive process. In what follows, we present a sequential algorithm 
to compute 0yK'''hich avoids matrix inversions. 

B. SEQUENTIAL LEAST-SQUARES ESTIMATION USING INSTRUMENTAL 
VARIABLES 

A sequential process for estimating the parameters of an unknown system requires 
fewer computations than a block process. In a manner similar to that presented in Hsia 
[Ref. 3: pp. 22-25] for the general least-squares case, the block IV estimation process 
described above can be converted into a sequential IV estimation process. Using the 
sequential process also allows the coeflicients to be updated based on the new data that 
becomes available. 

The derivation of the sequential estimation procedure consists of two parts. The 
first part is the derivation of an equation to update the data matrix, Q(/« + 1) , based 
on the previous data matrix, Q(/n) , and the new data: w{m ) , y{ni ) , and u{m -I- 1) 
where m represents the iteration. The second part involves developing an equation for 
updating the estimate of the parameters, 0,^{m+ 1), based on the previous estimate, 
On{m) , the previous data matrix, Q(/?i) , and the new data: w{m ) , y{m ) , and 

u{m -b 1) . 

Define the data matrix Q(/n) to be: 

Q(m) = [ZX]"’ (3-11) 

where Z„ is given by equation (3.5) and X„ is given by equation (2.8). The property of 
equation (3.3) assures that Q exists. Note that Q(/n) includes output data available 
through m and input data available through m + \ . Since both Z„ and X„ are 
m X (,V -b M) matrices, Q(wr) will be a (tV -b M) x (iV -b ^f) matrix. As the number of rows 
of Z and X increase to accommodate the increasing numbers of data points, the size of 
Q will remain the same. At the ne.xt sample time, i.e., at wi -b 1, the data matrix becomes: 

Q(m+l)-[zJ„,X„,.,r‘ (3.12) 

where the data matrices at -b 1 are given by: 
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(3.13) 



^wi+i 



z\m+ 1) 



V _ 
■^w+i ~ 



X„ 



x^(m + 1) 



(3.14) 



and i^{m+ 1) and x^(m+ 1) are vectors which contain the most recent data values. 
Substituting equations (3.13) and (3.14) into equation (3.12) and expanding, results in: 



Q(/n + 1) = 



[Z^ z(m + 1)] 



(3.15) 



xV+ 1) 



Expanding further yields: 

Q(/n + 1) = [Z^X„, + z(m + l)x^ + 1)]”’ (3.16) 

In equation (3.16) we see that two terms make up the new data matrix. The Z^„ term 
is all the data that was available through time m. The z{m + l)x^(m + 1) term contains 
the new data. To perform the inversion, let A = ZPC„, B = z(m+1), C=1 and 
D = \^{m + 1) . Then by the matrix inversion lemma: 

Q(/n+ 1) = A"’ -A"’B(C"' + DA"’B)"’DA"’ (3.17) 

Substituting the appropriate expressions for A, B, C, and D back into equation (3.17) 
yields the equation: 

Q(m + 1) = (Z^X„,)-’ - (ZXr‘z(^ + 1) 

. [l+xV+I)(Z^XJ-’z(m+l)]-’ (3.18) 

. xV + i)(zXr' 

Substituting Q(m) for (Z^XJ*' reduces equation (3.18) to: 

Q(/7j + 1) = Q(^n) — Q(ni)z(m + 1)[1 + x^(m + l)Q(m)z(/n + 1)]”’ 

. xV + DQM ' 
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This completes the first step of the derivation. Equation (3.19) expresses Q at time 
m + 1 in terms of the old Q and the new data. The term in the brackets is a scalar. 
Computational intensity has been reduced because a large matrix does not have to be 
generated and its inverse does not have to be calculated. 

Continuing with the derivation, the estimate 0;,. for data available through m can 
be written as: 



= (3.20) 

A 

The estimate d,y for data available through m + 1 can be written as: 

di,{m + 1 ) = (Z^+,X,+,)-'Z^+,y^+, (3.21) 

Substituting equation (3.12) into equation (3.21) results in an expression for the estimate 
of the parameters in terms of the new data matrix and all the available data given by: 

Bjiim + 1) = Q(m + l)Z^^iy„+, (3.22) 

5 m 

+ 1) = Q(m + l)[Z^ z(wi+l)] .... (3.23) 

_y{m + 1)_ 

+ 1) = Q(m + l)[Z^y^ + z(m + l)y(m + 1)] (3.24) 

Substituting for Q(m + 1) from equation (3.19) and expanding results in: 

Bilim + 1) = Q(^n)Z,^y^ 



- Q(m)z(OT + 1)[1 + x^(m + l)Q(m)z(/n + 1)] '.r^(m + l)Q(m)Z^y^ 

+ Q(m)z(m +!)>•(/«+ 1) (3.25) 

— Q(m)z(wi + 1)[1 + x^{m + l)Q(m)z(w + 1)]”' 

• \^{m + l)Q(OT)z(m + l)y(AM + 1) 

Although somewhat lengthy, this equation has the desired form. To simplify it, its last 
two terms can be arranged into the form: 

Q{m)z{m + 1)(1 — [1 + \^{m + l)Q(m)z(m + l)]~'x^(An + l)Q(m)z(A« + 1)} 

. y{m + 1) 
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Tlie terms within the braces can be thought of as the result of a previous application of 
the matrix inversion lemma with A~‘=l, B=l, C'' = 1 and 

D = x^(/n + l)Q(m)z(/n + 1 ) . Reversing the lemma results in: 

Q(m)z(m +!)[!+ x V + l)Q(/«)z(/?i + l)]~V("i + 1) (3.27) 

Replacing the last two terms in equation (3.25) with this result gives us: 

e,^m -r 1) = Q(m)Z^y^ - Q(/«)z(m +!)[!+ \^{tn + l)Q(m)z(m + 1)]"‘ 

. x\m + \)Q{m)Zly„ + Q(m)z(m + 1) (3.28) 

• [1 + x\m + l)Q(m)z(m + 1)]~V('” + 0 

Factoring Q(/^«)z(m + 1) and [1 + x^(/« + l)Q(m)z("?+ 1)]*' from the last two terms re- 
duces equation (3.28) further to: 

0,y{m -I- 1) = Qfm)Z^y„, -b Q(m)z(m -b 1)[1 -b x^{m -b l)Q(m)z(w -b 1)]"’ ^ 9 ^ 

. [y{m -b I) - X V + l)Q('”)Z^y J 



Substituting equation (3.1 1 ) into equation (3.20) and then equation (3.20) into equation 
(3.29) yields the final form for the update of the estimate of the parameters: 



-b 1) = Oiiim) + Q(w)z(/n -b 1) 

• [1 -b X^(/?I -b l)Q(/?l)z(/7I -b !)]”'[>■(/« -b 1) — xJ{m -b l)0/t/(«l)] 



(3.30) 



This is the desired result for updating the estimate of the parameters. Note that like 
equation (3.19), the matrix inversion of equation (3.21) has been reduced to inversion 
of a scalar. Equation (3.30) describes the update of d,^m+ 1 ) in terms of the previous 

A 

estimate of the parameters, Qnim), the previous data matrix, Q(m), and the new data: 
vv(/m) , y{>n) , and u[m + 1 ) . 

C. TESTING THE SEQUENTIAL INSTRUMENTAL VARIABLE ALGORITHM 
Equations (3.19) and (3.30) above comprise the sequential IV algorithm. Several 
tests of this algorithm were made using second and third-order filters as unknown sys- 
tems. Tests were run via computer simulation using the filters to generate the output 
data. A Gau.ssian random process with zero mean and unit variance was used as the 
input. The input was produced by IMSL subroutine GGNML. Graphs were created 
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using DISSPLA. Table 1 on page 17 shows pole and zero locations as well as numera- 
tor and denominator parameters for the test filters. 



Table 1. TEST SYSTEMS FOR THE IV MODETING METHOD 



TEST 


LOCATIONS 


LOCATIONS 


AR 


.MA 


FILTER 


OF POLES 


OF ZEROS 


PAR.AMETERS 


PARAMETERS 


T2 


0.445 + j0.228 


0.4 + jl.273 


1.0 


0.5 




0.445 - jO.228 


0.4 - jl.273 


-0.89 


-0.4 




0.25 


0.89 


T2X 


0.445 + j0.228 


0.4-t-jO.S 


1.0 


1.0 




0.445 -jO.228 


0.4 -jO.8 


-0.89 


-0.80 




0.25 


0.80 


T3 


0.6605 


-1.0 


1.0 


0.0154 




0.6647 -)- j0.o02 


-1.0 


-1.99 


0.0462 




0.6647-j0.502 


-1.0 


1.57 


0.0462 






-0.458 


0.0154 



Results of the tests are sho\sm in graphical form in Figure 4 on page 18, Figure 5 
on page 19. and Figure 6 on page 20. Dashed lines indicate the true values of the pa- 
rameters. Solid lines are the IV method's estimates. 

For both second-order test cases shown, the algorithm converged quickly and 
produced accurate results. For the third-order test case, convergence took longer but 
the values were accurate. A third-order system is more comple.x than a second-order 
system, so conceivably it would require more iterations to converge. The number of it- 
erations required is of the same order as the method of ordinary least-squares. 

Table 2 on page 21 contains the IV algorithm's best estimates of the parameters and 
the number of iterations required to converge to those estimates. It also shows the ab- 
solute and percent errors from the true parameters. 
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(A) 




(B) 



Figure 4. Second-order test case T2. (A) MA parameters. (B) AR parameters. 
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PARAMETER VALUES PARAMETER VALUES 



1.5 




6 2 4 6 0 lb 12 14 16 lb 20 



ITERATIONS 



(A) 




ITERATIONS 



(B) 

Figure 5. Second-order test case T2N. (A) MA parameters. (B) AR parameters. 
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PARAMETER VALUES 




ITERATIONS 



(A) 




(B) 



Figure 6. Third-order test case T3. (A) MA parameters. (B) AR parameters. 
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Table 2. COEFFICIENT ESTIMATES BY THE IV MODELING METHOD. 



TEST 

FILTER 


PAR.AMETER 

ESTIMATE 


ABSOLUTE 

ERROR 


PERCENT 

ERROR 


ITEIUATIONS 


T2 


0.500 


0.0 


0.0 


10 




-0.396 


-b 0.004 


O.IO 






O.SSS 


-0.002 


■ 0.22 






1.000 


0.0 


0.0 






-O.SSS 


+ 0.002 


0.22 






0.244 


-0.006 


2.40 




T2N 


1.000 


0.0 


0.0 


10 




-0.794 


+ 0.006 


0.750 






0.794 


-0.006 


0.750 






1.000 


0.0 


0.0 






-0.SS7 


+ 0.003 


0.34 






0.243 


-0.007 


2.80 




T3 


0.0154 


0.0 


0.0 


1000 




0.0466 


+ 0.0004 


0.87 






0.0475 


+ 0.0013 


2.S1 






0.0169 


+ 0.0015 


9.74 






1.000 


0.0 


0.0 






-1.^)6 


-0.03 


3.0 






1.532 


-0.040 


2.01 






0.4379 


-0.0204 


4.45 





IV. SYSTEM IDENTIFICATION USING AN ITERATIVE 
MULTICHANNEL APPROACH 



A. INTRODUCTION 

This chapter presents an alternate system identification method, the iterative multi- 
channel approach. This approach differs from the IV method and the method of ordi- 
nary least-squares presented in the preceding chapters in that it processes the input and 
output data from the unknown system in separate channels. In its block processing 
form one advantage over the IV and ordinary least-squares methods is a reduction in the 
sizes of the data matrices. As a result, the computational complexity of the multichannel 
algorithm is on the order of -I- , where .Vf is the order of the .\IA part and A' is the 

order of the AR part. In contrast, the block IV and ordinaiy' least-squares methods re- 
quire computations on the order of (.1/ + i\y . 

B. PREVIOUS MULTICHANNEL METHODS 

Whittle [Ref 6: pp. 129-130] was the first to develop a multichannel solution for the 
.ARMA modeling problem. He sought to extend the recursive Durbin solution for esti- 
mating the parameters of a single variable autoregressive process to a multivariable 
autoregressive process. He discovered that to do this he would have to fit the data to 
two autoregressive processes simultaneously. One of the autoregressions would use 
present data samples to predict the value of the data one time step in the future. This 
is called forward prediction. The second autoregression would use present data samples 
to predict the value of the data at the previous time instant and is referred to as back- 
ward prediction. Sometime during this research. Whittle determined that if the input 
was derived from a MA scheme, making the process ARMA, then the solution would 
remain the same provided the correlations of the input used in the parameter estimation 
process had shifts greater than the MA scheme. Whittle's use of the two separate and 
simultaneous autoregressions to model an ARMA process can be thought of as a 
multichannel modeling approach. 

Further work in the area of multichannel ARMA modeling was conducted by Perr>' 
and Parker [Ref 7; pp. 509-510]. They started out with the ARMA problem formulation 
discussed in Chapter 2. Using the method of ordinarx' least-squares to minimize the 
mean square error, they found the solution for the estimate of the parameters to be the 
Wiener solution given by: 
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( 4 . 1 ) 



1 




■ b' 




1 

1 


cc: 
— 1 


K'w, 


_a'_ 




1 



In equation (4.1) is a matri.x of autocorrelations of the past outputs. R„.„, is a matri.x 
of autocorrelations of the inputs, R^„. and R„.^ are crosscorrelations of the input and 
output data, b is a vector of AR parameters, and a' is a vector of MA parameters. In 
addition, is a vector of autocorrelations of past output data with the current output 
and r^„. is a vector of crosscorrelations of input data with the current output. By as- 
suming the first .MA parameter, a\ was known, they were able to treat it as a gain and 
factor it out of all the other MA parameters. This allowed them [Ref. 7: pp. 509-510] 
to extract the (:V-t- 1)” row and column of equation (4.1) and rewrite the solution in the 
form: 



1 


^yu 


b' 




~Vyy' 




'g/ 


1 


Ku_ 


a_ 








Juu. 



In equation (4.2), a' has been rewritten as a to indicate that the a\ term has been fac- 
tored out of all of the MA parameters. Reasoning that equation (4.2), the AR.MA sol- 
ution, was a generalization of the AR solution, they figured that it must have a recursive 
solution consisting of some combination of the Levinson-Durbin algorithm, a recursive 
solution for the AR problem. They then determined equation (4.2) was in a form similar 
to Whittle's formulation of the problem. So they reasoned that they could use a form 
of Whittle's solution to solve the .A.RM.A modeling problem. Like Whittle, their solution 
consists of a forward and a backward autoregression. It uses two coupled lattice filters 
to process the input and output data. Off diagonal elements of the lattice coefficient 
matrices specify the coupling points of the two lattices. 

C. ITERATIVE APPROACH TO MULTICHANNEL ARMA MODELING 

This thesis proposes another solution to the AR.MA modeling problem using the 
multichannel approach. It is an iterative approach with no direct coupling of the two 
channels. However, note that there is an implicit coupling in the sense that the AR.M.A 
system output samples >’(/?) are a function of the present and past input samples u{n), 
and the past output samples. This is shown in equation (2.1). The approach proposed 
uses two separate finite-impulse-response (FIR) filters to estimate the unknown system. 
One filter estimates the AR part of the unknown system. The second estimates the .MA 
part of the unknown system. Figure 7 on page 25 shows the structure of this approach. 



23 



The derivation of the equations for this approach follows the method of ordinary' least- 
squares. 

A(z) 

From Figure 7 on page 25, the value of the signal Y{z) is seen to be ^ • 

When this signal passes through C(z) , if C(r) is close to B(r), the resulting signal X{z) is 
approximately A(z)(/(z). Also from Figure 7 on page 25, the value of Z(z) is seen to be 
A(z)C*(2) provided D(z) is close to A(z). The difference of these two signals forms the 
error which we seek to minimize by the method of ordinary least-squares. In minimizing 
the error we seek to drive D(z) and C(z) as close as possible to A(z) and B(z) , respec- 
tively. 

Referring to Figure 7 on page 25, signals z{n) and jr(«) are defined as the outputs 



from two FIR filters and are given by the equations: 

z(w) = dQu{n) + d^u{n — 1) + d 2 U{n — 2J -I — -I- dj^ju{n — ^^) = u^(«)d (4.3) 

=T(«) + <^i>'(" - 1) + ^’xv(rt - 2) -I- ••• -I- c^'{n - N) = y^(«)c (4.4) 

where the vectors d and c represent the systems D(z) and C(z), respectively. The vectors 
d, u(«), c, andy(/i) are given by the following equations: 

d = [4 d, d2 ... (4.5) 

u(«) = [u(/i) !/(/;— 1) u{n — 2) ... «(« — i\/)]^ (4.6) 

c = [l c, cj ... c,v]^ (4.7) 

y(A?) = [>•(«) y(/i - 1) y(« - 2) ...y(i; - (4.8) 



The d parameters are estimates of the MA portion of the ARMA process. The c pa- 
rameters approximate its AR portion. The vector u(/?) is the input data vector of length 
M, the order of the M A part, and y(«) is the output data vector equal of length N, the 
order of the AR part. 

Forming the error between x and z results in: 

e{n) = z(«) — jf(«) = u^(«)d — y\n)c (4.9) 

The least-squares performance criterion is: 
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UNKNOWN SYSTEM 




e(n) = z(n) - x(n) 



Figure 7. Multichannel modeling process 



L L 

J = ~ 

/7=1 /2=l 



(4.10) 



Substituting for e[n) from equation (4.9), we can write equation (4.10) in vector form as: 



L 

y=^(u^d-yV 

/ 1=1 



(4.11) 



where we have dropped the time index, n, for convenience. Expanding equation (4.11) 
results in: 



L 

J = ^d^uu^d + c^yj’^c — 2d^uy^c 



(4.12) 



n=\ 
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We notice that the performance criterion is a function of both d and c. Minimizing the 
performance criterion by differentiating with respect to the vector c and equating the 
results to zero yields; 



L L 

-2^(yu^)d (4.13) 

/l=l .«=! 

Similarly, differentiating the performance criterion with respect to the vector d and 
equating the result to zero yields: 



dJ 

PC 



= 0 = 0+2 



Vot 



)c 



iL 

cd 



= 0 



L 



L 




(uu^)d - 2 




(uy^)c 



? 2 =\ 



n=l 



(4.14) 



Solving equation (4.13) for c and equation (4.14) for d results in two equations for 
estimating the AR and MA parameters of the unknown system given by: 



c 



= R^„d 



and 



d = RjRu>c 



where 



Ruu = ^U«^ = 

rt=l 



^uuO) >uui^) 






''uuU’^’) 



■ I'uM 



(4.15) 



(4.16) 



(4.17) 



is an estimate of the input autocorrelation matrix. The elements of R^„ are computed 
using the unbiased method as follows: 



L-l 

„(0 = - 0 for / = 0, 1, 2, ... , 

j=0 



M 



(4.18) 



26 



Matrices R„^, and R^„ appearing in equations (4.15) and {4. 1 6) have structures iden- 
tical to equation (4.17), where R^.^ is the estimate of the output autocorrelation matrix, 
and R^j and R^„ are estimates of the crosscorrelation matrices. Note that R^,,, = R^^. . The 
elements of these matrices; and are computed as follows; 

L-l 

- /) for / = 0, 1 , 2, ... , iV (4. 1 9) 

MO 

L-l 

^ for ^ = 0. 1. 2, ... ,M (4.20) 

y=o 



and 



L-l 



VW = y^^‘(/>(/‘ - 0 for / = 0, 1, 2, 

MO 



N 



(4.21) 



Up to this point, following the standard least-squares procedure has resulted in two 
dependent or coupled equations to solve for the parameters of an unknown system 

modeled as an ARMA process. How best to solve these equations? By iteration. The 

steps of the iterative process are to 

• Calculate the correlation matrices and vectors from the available data. 

• Initialize c . Exploit the fact that the first parameter in c is, Cg = 1. 

• Solve for d from this initial c . 

• Solve for c from d . 

• Repeat beginning at the third step. 



Here is a summary of the equations in proper order for implementing the algorithm: 

• Compute R„„ , R^^ and R„^ from equations (4.17) to (4.21). Note that R,.„ = 

• Initialization; 

c^°^ = [l 0 ... 0]^ (4.22) 



• For k = 0lo K 



d^"-^'> = R-'R„/"> 



(4.23) 

(4.24) 



where k is the iteration number. 
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This is a ven’ simple algorithm. For the case where the AR and MA orders are equal, 
the correlation matrices are half the size of the block data matrices which must be gen- 
erated and inverted in the IV algorithm. 

1. Testing the Multichannel Iterative Algorithm 

We tested the algorithm by computer simulation using second and third-order 
filters as unknown systems. Table 3 shows pole and zero locations as well as MA and 
.AR parameters for the test filters. Data lengths of 500 data points were used to calculate 
the autocorrelation and crosscorrelation matrices. Besides the reported cases, we tested 
the algorithm on several other second, third and nii.xed-order cases. The performance 
of the algorithm in all cases that we tested was fairly uniform. 



Table 3. TEST SYSTEMS FOR ITERATIVE MULTICHANNEL MODELING 
METHOD 



TEST 


LOCATIO.X 


LOCATION 


AR 


MA 


FILTER 


OF POI.ES 


OF ZEROS 


PA R,A METERS 


PAICAMETERS 


T2 


0.445 + j0.228 


0.4 + jl.273 


1.0 


0.5 




0.445 - j0.22S 


0.4 -j 1.273 


-0.89 


-0.4 




0.25 


0.89 


T2X 


0.445 A jO. 228 


0.4Aj0.8 


I.O 


1.0 




0.445 -jO.228 


0.4 -jO.S 


-0.89 


-0.80 




0.25 


0.80 


T3 


0.6605 


-1.0 


1.0 


0.0154 




0.6647 -^j0.502n 


-1.0 


-1.99 


0.0462 




0.6647 -jO.5020 


-1.0 


1.572 


0.0462 






-0.4583 


0.0154 



Results of the tests are shown in graphical form in Figure 8 on page 29, 
Figure 9 on page 30. and Figure 10 on page 31. Dashed lines indicate the true values 
of the parameters. Solid lines show the values the algorithm calculated. 

Table 4 on page 32 contains the algorithm's best estimates of the parameters, 
along with the number of iterations required to converge to those estimates. It also 
shows the absolute and percent errors from the true parameters. For the second-order 
test cases, convergence to the true parameter values occurred within 20 iterations. The 
third-order test case took 14 iterations to converge to its steady-state values; however, 

these values were not the true parameter values. 
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PARAMETER VALUES PARAMETER VALUES 



1.0 




ITERATIONS 



(A) 




(B) 



Figure 8. Second-order test case T2. (A) MA parameters. (B) AR parameters. 
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PARAMETER VALUES PARAMETER VALUES 




(A) 




(B) 



Figure 9. Second-order test case T2N, (A) MA parameters. (B) AR parameters. 
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PARAMCTRR VAI.UES PARAMETER VALUES 



0.3 




- 0.1 ^ ^ ^ : : : ^ : ' 

0 2 4 6 8 10 12 14 16 18 20 

ITERATIONS 



(A) 




(B) 



Figure 10. Third-order te.st case T3. (A) MA parameters. (B) .AR parameters. 
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Table 4. PARAMETER ESTIMATES BY THE ITERATIVE MULTICHANNEL 
ALGORITHM. 



TEST 

FILTER 


PA R.A METER 
ESTIMATE 


ABSOLUTE 

ERROR 


PERCENT 

ERROR 


ITER^ATIONS 


T2 


0.500 


0.0 


0.0 


5 




-0.393 


-0.007 


1.75 






0.890 


0.0 


0.0 






LOGO 


0.0 


0.0 






-0.SS9 


+ 0.001 


0.11 






0.247 


-0.003 


1.20 




T2N 


LOGO 


0.0 


0.0 


20 




-0.794 


+ 0.006 


0.75 






0.79S 


-0.002 


0.25 






1.000 


0.0 


0.0 






-0.886 


+ 0.004 


0.45 






0.245 


-0.005 


2.00 




T3 


0.0153 


-0.0001 


0.65 


14 




0.0487 


+ 0.0025 


5.41 






0.0590 


+ 0.0128 


27.71 






0.0287 


+ 0.01 33 


86.36 






0.99 


-0.01 


1.0 






-1.79 


+ 0.20 


10.05 






1.267 


-0.305 


19.40 






-0.3O86 


+ 0.1497 


32.66 





2. Stopping Parameter 

In tests of third-order systems, the parameter estimates swung through or close 
to the true coefficient values and converged to values somewhat removed from the true 
values. We developed a stopping parameter to Hag the iteration when the estimates were 
closest to the true values. This occurs when the error is smallest. Referring to 
Figure 7 on page 25, if D(r) is equal to A(x) and C(c) is equal to B(r), .y and z will both 
equal A(z)L'(-). The error will be zero. The farther removed D(z) and Cfz) are from 
A(z) and B(z), respectively, the larger the error becomes. 

The stopping parameter is calculated from the difference of the z and x values 
at every' iteration. After the parameter vectors d and c are estimated for a particular it- 
eration, the stopping parameter is calculated from: 

^ki") = - ^ki") = yl^^k - ^k'^k 



32 



where k is the iteration number and d . u, c and y are defined in equations (4.5) to (4.S). 
The parameter vectors d and c represent the systems D(c) and C(z), respectively. 

Figure 11 on page 34 and Figure 12 on page 35 graph the stopping parameter 
(dotted line) along with the estimated and true values of the parameters. They show that 
when the stopping parameter is smallest, the parameters are closest to their true values. 
Table 5 shows the improvement in the estimates of the parameters resulting from 
choosing those values when the stopping parameter is smallest. 



Table 5. PARAMETER ESTIMATES WHEN STOPPING PARAMETER IS 
SMALLEST, 



TEST 

FILTER 


PAR.AMETER 

ESTl.MATE 


ABSOLUTE 

ERROR 


PERCENT 

ERROR 


ITER.ATIONS 


T3 


0.0156 


+ 0.0002 


1.30 


10 




0.0478 


+ 0.0016 


3.46 






0.0536 


+ 0.0074 


16.02 






0.0196 


+ 0.0042 


27.27 






0.97 


-0.03 


3.0 






-1.84 


+ 0.15 


7.54 






1.375 


-0.197 


12.53 






-0.3672 


+ 0.0911 


19.88 





The stopping parameter can be used in a real modeling situation because it 
comes from the data and the estimates of the parameters. Another measure of how well 
the estimates of the parameters fit the actual system is the norm of the coefficient error. 
This cannot be used in a real modeling situation however, because the values of the true 
parameters are not known. W’e calculated it for the test cases as a check on the appro- 
priateness of using the stopping parameter. Figure 13 on page 36 and Figure 14 on 
page 37 graph the stopping parameter (dotted line) and the norm of the coefficient error 
for test cases T2 and T3. On both graphs the two curves correspond well. Both reach 
their minimum value at the same point, the point where the estimates of the parameters 
are closest to their true values. 

3. Linear-prediction of the Denominator Coefficients 

The iterative algorithm detailed in equations (4.22) to (4.24) starts by initializing 
the .AR parameter estimates to: 

c^^> = [l 0 ... 0]^ (4.26) 
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PARAMETER VALUES PARAMETER VALUES 



1.0 



0.15 




(A) 




ITERATIONS 



(B) 

Figure 11. Stopping parameter example for test case T2. 



0.10 



0.05 



0.00 
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STOPPING PARAMETER STOPPING PARAMETER 



PARAMETER VALUES PARAMETER VALUES 




(A) 




ITERATIONS 



(B) 



Figure 12. Stopping parameter example for test case T3. 
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STOPPING PARAMETER 




Figure 13. Stopping parameter and norm of the coefficient error for test case T2. 

where Cq is known to be 1 from the Z transform of the ARM A difference equation, 
equation (2.1). The Z transform is given by; 



] (zj[l -1- C,Z ’ -t- C22 


+ •••]= + - 3 


(4.27) 


ri--) 


t/fj -(- d^z 


(4.28) 


l\z) 


1 + CiZ + C 2 Z + "• 



The initial estimates for the other AR parameters are zero. This can be far from 
their actual values. A closer estimate of the other AR parameters should result in 
quicker convergence for all parameters. A closer estimate of the AR parameters can be 
obtained by using linear-prediction techniques. Figure 15 on page 39 shows the ap- 
proach used. In Figure 15 on page 39, y(«) is the output from the unknown system. 
The system C'(r). which is represented by the vector c' , is the linear-prediction filter used 
to estimate y(//). It uses the previous n — .V samples of the output to generate a current 
estimate. This is given by the equation: 
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Figure 14. Stopping parameter and norm of the coefficient error for test case T3. 

j'-(rt) = y^(n- 1)C' (4.29) 

where c' is a vector of the tap weights of the autoregressive process given by: 

c' = [c'o c', ...c'„_.nT (4-30) 

and y(/7 — 1) is a vector of the past output values given by: 

y(// - 1 ) = iv(n - 1) y{n - 2) ... >•(« - .\-)]^ (4.31) 

Following least-squares techniques, we form the error between the estimate and the ac- 
tual value of the output. The sum of the squares of the errors becomes the performance 
criterion. This is differentiated with respect to the tap weights and set equal to zero. 
Solving this for the tap weights results in the equation: 
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1 his is the standard Weiner filter solution [Ref 8; p. 32]. It tells us that the best estimate 
of the AR parameters can be found from the correlations of the output data. The matrix 
is the autocorrelation matrix of the past outputs and is autocorrelation vector of 
the past outputs with the current output. In all cases tested, we did not achieve any 
significant improvement in the accuracy of the estimates of the parameters, or in the 
speed of convergence, using the straight linear prediction of equation (4.32). 

A modification to this approach, which we refer to as modified linear-prediction, 
uses correlation lags beginning on the order of the MA portion of the ARMA process. 
For example, correlations for calculating for a third-order system would start at a lag 
of three and increase to a lag of five. Correlations for calculating would start at a lag 
of four and increase to a lag of six. This ensures that the effect of the MA part of the 
unknown system is removed from the linear-prediction of the AR part. This modified 
method of linear-prediction is given by the equation: 



c'o 

c'l 


= 


fyyiq) 

^yyi.R ^yyi^) 


• ryy(^-p+2) 


ryy{q + 1) 
ryyiq -1- 2) 


(4.33) 


c'„-,v 




>'yyi(} + P-^ fyyiq + p-2) . 


ryyiq) 


fyyiq + p) 





where q is the order of the MA portion and p is the order of the /\R portion. [Ref 2: 

p. 182] 

Figure 16 on page 40 shows the results of using modified linear-prediction with 
third-order test case T3. When comparing this graph to the estimates obtained without 
linear-prediction, shown in Figure 12 on page 35, note that the vertical axes have dif- 
ferent scales. Table 6 on page 39 lists the values of the estimates at iteration 10 and 
compares them with the true values via the absolute and percent errors. A comparison 
of Table 6 on page 39 with Table 5 on page 33, the best estimates without the use of 
modified linear-prediction, shows that modified linear prediction has significantly in- 
creased the accuracy of the AR estimates at the tenth iteration. The accuracy of the 
MA estimates remains approximately the same. The tenth iteration was chosen as the 
point to select the parameter values because in both cases this was the iteration where 
the stopping parameter had the smallest value. 
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y(n) 




Figure 15. Linear prediction block diagram 



Table 6. PARAMETER ESTIMATES USING MODIFIED 

LINEAR-PREDICTION FOR INITIAL ESTIMATE OF AR PARAME- 
TERS. 



TEST 

FILTER 


PA R.A METER 
ESTIMATE 


ABSOLLIE 

ERROR 


PERCENT 

ERROR 


1TER.AT10NS 


T3 


0.015(3 


-r 0.0002 


1.30 


10 




0.0-!S5 


-^0.(j023 


4.98 






0.')5 1 2 


0.0050 


10.80 






0.0101 


-i- 0.0053 


34.42 






l.OO 


0.0 


0.0 






-1.98 


+ 0.01 


0.50 






1.553 


-0.019 


1.21 






-0.4458 


+ 0.0125 


2.73 





D. FORMULATION OF THE SEQUENTIAL MULTICHANNEL APPROACH 
To decrease the computational intensity of updating the estimates of the .AR and 
-M.A parameters due to new data, an algorithm to sequentially process the data has been 
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Figure 16. Parameter estimates for test case T3 using modified linear-prediction for 
initial estimate of AR parameters. 
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developed. Development begins with the performance criterion seen previously for the 
block data case in equation (4.10). From that starting point, equations are developed 
which relate new estimates of the parameters to the previous estimates and the new data. 
Separate but similar equations are developed for the MA and the AR coefficients. Due 
to the similar nature of the development of these equations, only the development of the 
equations for the AR coefficients is presented here. However, the final results for both 
the AR and MA coefficients are given. 

The performance criterion can be written as: 

n ^ /I 



E.xpanding this results in: 



J = ~ (435) 

i=0 

where c and y are defined in equations (4.7) and (4.8) and z is the scalar signal at the 
output of D. Differentiating the performance criterion with respect to c and setting the 
result equal to zero yields: 



i=G i=0 



Equation (4.36) can also be written as 



n n 



Yj^iyi= 

1=0 i=0 



Solving for the AR parameter vector results in: 



n n 



/=o /=o 



(4.36) 



(4.37) 



(4.38) 
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Because c is an estimate based on data available through time n we signily this by in- 
troducing the index « on c to yield: 



n n 

c, = (Z,^-,yIr''Z^,y, ( 4 . 39 ) 

/=0 i =0 

The first step in formulating the sequential algorithm is to develop an update 
equation for the estimate of the AR parameters. Applying the method presented in 
Graupe [Ref 9: p. 124], vve first define a new matrix P~‘ as 

p;' - fj-,yl (4.40) 

(=0 

This is a matrix of the output data of the unknown system. At the previous time n — 1 
this matrix is written as: 



p;-i = Ey/jf (4.41) 

/=0 

By substituting equation (4.40) into equation (4.39) the estimate of the AR parameters 
can be rewritten as: 



C, - P„X!"<yi (■'.42) 

/■=o 

The right side of equation (4.42) needs to be converted into an expression containing the 
previous estimate of the parameters plus a correction term. It needs the past value of 
c„ which is c„_j. From equation (4.39) we can write: 

n-l n-\ 

C„-, = (Ewri-'Ewi (4.43) 

/=0 /=0 

This can be rewritten as: 

^-1 n-\ 

= iy^yiy[)Cn-\ (4-44) 

2=0 /=0 
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Premultiplying equation (4.42) by and then separating the last term from the sum- 
mation results in: 



n-\ 

/=0 

r?-l 

Substituting for X^.y, in equation (4.45) from its equivalent expression in equation (4.44) 

1=0 

yields: 



n-l 

’c = (^y,y/)c„_, + 2„y„ (4.46) 

i=0 

By adding y„yfc„_, to the right side of equation (4.46) and grouping it with the summa- 
tion, the summation will range from / = 0 to « . In order not to affect the value on the 
right side of equation (4.46), y„y„^c„_, must also be subtracted from the right-hand side, 
which yields: 



n-l 



(/ )^n— 1 "b "b ^ nb n^n—\ 1 



i=0 



Combining y„yjc„_, with the summation as describe above results in: 



(4.47) 



1^ /3 i 1 "b ynJfi^n— 1 ("^'^S) 

1=0 

n 

Replacing Xyjf witli its equivalent expression from equation (4.40) yields: 
t=r 

'c„ = P« 'c„_i + y„(z„ - y„^c„_i) (4.49) 

Premultiplying by P„ results in the following equation for the update of the estimate of 
the AR parameters: 

c„ = c„_i -f P„y„(z„ - y^:„_i) (4.50) 

This is the desired result. It relates the estimate of the parameters at time A' to the es- 
timate at the previous time. A' — 1, plus the new output data vector, y„, and the error at 
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lime .V. Error is represented by the term. The corresponding equation for the 

MA parameters is; 

+ Q,jU,j(vV„ — u„d^_]) (4.51) 

In equation (4.51), Q is a matrix of the input data of the unknown system given by: 

n 

Q;' = ^u,uf (4.52) 

i=0 

Finally, we need a sequential update for P„ and Q„ to complete the sequential algo- 
rithm. This is accomplished by using a form of the matrix inversion lemma. 

By pulling the last term out of the summation, the definition of given in equation 
(4.40) can be rewritten as; 



P«‘ = ^yi>’r+ynyr 

i=0 



(4.53) 



n-I 



Substituting for Xy/jT its equivalent expression from equation (4.41) results in 

/=0 

P“' = p“’ J. Y 

* /2 * ;?-l ^ Jnin 



(4.54) 



Inverting both sides of the equation results in: 



= (p;l, + y,, >■!)■' 



(4.55) 



Let A = P;i|, B = y„, C = 1 , and D = y^. Then, by the matrix inversion lemma; 

P„ = - A“’B(C~’ + DA“’B)“’DA"' 



(4.56) 



Substituting the appropriate expressions into equation (4.56) results in: 

p, - (p;ii)'' - (p;-,r'y,[i + )f(p;iir‘yj'‘yr(p;'i)' 



(4.57) 



This reduces to; 



p„ = p„_, - P„_,y„[i + y«P„-iyJ 'y,[p«_i 



(4.58) 



Using this same procedure, the update for Q„ is; 
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(4.59) 



Q„ = Q„_, - Q„_iU„[l + uJ'q„_,uJ ‘uJq„_, 

A reduction in the computational intensity has been achieved by reducing the matrix 
inversion of equation (4.55) to inversion of a scalar in equation (4.57). Inversion of 
these scalars is much simpler than inversion of the Rj,^ and matrices of the block 
processing case. 

The sequential multichannel algorithm is summarized below: 



• The parameter update equations: 




~ ^n—\ ~ y/7^/7— l) 


(4.60) 




(4.61) 


• The update equations for the P and Q matrices: 




P /7 “ P/7-1 ““ P/ 7 -iy/X^ T y /7 P/7-1 )/?] y /7 P/7-1 


(4.62) 


Q„ = Q„_, - Q„_iU„[l + 


(4.63) 



The reduction in computational intensity comes with a trade-off. Now the algo- 
rithm is more complex to use. Updates are required for P and Q as well as c and d where 
before, in the block multichannel algorithm, updates were only required for c and d. But. 
as in the sequential IV case, an added advantage of the sequential multichannel algo- 
rithm is it allows updates of the estimates of the parameters based upon new data. 
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V. SUMMARY 



In this thesis we set out to develop two algorithms for modeling unknown systems 
as ARMA processes. These are the IV method of system identification presented in 
Chapter 3, which is a modification of the method of ordinary least-squares, and the it- 
erative multichannel method presented in Chapter 4. 

The IV method is not a new concept in either its block or sequential processing 
forms. However, our derivation of the sequential algorithm was done independently of 
other IV sequential algorithms. We chose the IV method because it reportedly has good 
noise handling capabilities and yields consistent and unbiased estimates of the unknown 
system’s parameters. It also remains as easy to use as the method of ordinary least- 
squares. W’e also wanted to gain familiarity with it because it was a possible candidate 
for incorporation into the multichannel method. 

Operating alone, the IV method produces accurate estimates of the unknown sys- 
tem's parameters. Convergence was within 20 iterations for several second-order sys- 
tems that we tested. Convergence slows down as the system order increases. However, 
the results do converge to the actual system parameters given sufficient processing time. 
The performance of the IV method is similar to the performance of the method of ordi- 
nary least-squares. 

The proposed iterative multichannel algorithm is new in both its block and sequen- 
tial processing forms to the best of our knowledge. It is ver>‘ simple to use in the block 
form. It achieves accurate results for second-order systems but worse results for third- 
order systems with block correlation elements calculated based on only 500 data points. 
Implementing the stopping parameter increases the accuracy when the parameter esti- 
mates converge but not to the true parameters. Due to its ability to separately process 
the input and output data from the unknown system, correlation matrices in the multi- 
channel block processing case are half the size of correlation matrices required for the 
single channel block processing case. This feature reduces the computational intensity 
over what is required for the conventional least-squares processing case. The number 
of iterations required for convergence seems to be independent of the order of the sys- 
tem. However, the accuracy of the estimates suffer as the order of the system increases. 
Using linear prediction to estimate the initial values of the AR parameters did not speed 
up convergence or increase the accuracy of the parameter estimates. However, using 
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modified linear prediction significantly increased the accuracy of the AR parameter es- 
timates. although it had no effect on the MA parameter estimates. 

We formulated the multichannel sequential algorithm. This allows the estimates of 
the parameters to be updated as new data becomes available. But we have not tested 
this algorithm. It needs checking using a variety of second and third-order test cases to 
verify that it works. During testing, guidelines need to be developed for the best meth- 
ods to initialize the P and Q matrices to achieve the quickest convergence and the most 
accurate parameter estimates. 

As mentioned above, one of the reasons for investigating the IV method of system 
identification was for possible inclusion into the multichannel algorithm, the hope being 
that the favorable noise performance of the IV method would improve the performance 
of the multichannel method. This is another area that remains unexplored. 

The block multichannel and IV methods achieved similar results for second-order 
test cases. Convergence to the actual system parameters came within 20 iterations for 
both algorithms. However, for third-order systems, convergence was much quicker with 
the multichannel block algorithm than with the sequential IV algorithm. But the pa- 
rameter estimates by the IV method were more accurate than by the multichannel block 
method. A combination of the two algorithms has the potential for incorporating their 
unique advantages into a better overall parameter estimation method. 

Areas for further research are listed below; 

• Verify that the multichannel sequential algorithm developed in Chapter 4 works as 
a means of modeling an unknown system as an AR.VIA process. 

• Investigate the possibility of incorporating the IV method into the multichannel 
sequential algorithm. 

• Analyze why initializing the AR parameters to the values calculated by linear pre- 
diction improves the speed of convergence of the AR parameters in the multi- 
channel block algorithm but does not improve the convergence of the .MA 
parameters. Identify a method for obtaining an initial estimate of the .MA param- 
eters to improve their speed of convergence and accuracy. 

• Investigate the effects of increasing the number of data points used to calculate the 
correlation matrices for the multichannel block algorithm on the accuracy of the 
parameter estimates and their speed of convergence. 

• Investigate the performance of the IV method with noise present. Compare this 
performance with the performance of the method of ordinary least-squares with 
noise present. 
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APPENDIX 



A. INSTRUMENTAL VARIABLE ALGORITHM 



* * 



* PAUL DAL SANTO * 

* IV ALGORITHM * 

* 4/12/88 * 

^ Vr 

* THIS PROGRAM CALCULATES THE AR AND MA PARAMETERS OF A * 

* TEST SYSTEM BASED UPON ITS INPUT AND OUTPUT DATA * 

* BY USING THE SEQUENTIAL IV METHOD. * 

Vf * 

* * 



, 



■jVyrVf*?Wf*yf^V'sV**Vc’**Vf^Wf*yf*yf^VVfVf^V*V«’**?V’sVVfyr*':^ <**■;??* 






VARIABLE DEFINITIONS **** 



* 

it 



it 

y. 

*V 



* 

it 



it 

it 

it 



it 



RIVCOF ARRAY FOR STORING THE AR AND MA PARAMETERS 
THE PROGRAM CALCULATES 

Z VECTOR OF DATA FROM THE OUTPUT OF THE AUXILIARY 

MODEL AND THE INPUT TO THE TEST SYSTEM 
Z TPO TRANSPOSE OF VECTOR Z 

X VECTOR OF DATA FROM THE OUTPUT AND INPUT OF THE 

TEST SYSTEM 

X TPO TRANSPOSE OF THE X VECTOR 

U STORAGE FOR INPUT DATA 

Y STORAGE FOR OUTPUT OF TEST SYSTEM 

W STORAGE FOR OUTPUT OF THE AUXILIARY MODEL 

QMAT THE Q MATRIX OF THE IV ALGORITHM 

IV VECTOR OF PARAMETERS CALCULATED BY THE ALGORITHM 

COR RESULT OF INTERMEDIATE STEP IN ALGORITHM CALCULATION 

COEF VECTOR OF TRUE PARAMETERS OF TEST SYSTEM 

SCALAR RESULT OF SCALAR INVERSION IN INTERMEDIATE STEP COR 
SCALR2 INTERMEDIATE STEP WHEN CALCULATING THE NEW IV VECTOR 
SEED INITIALIZATION FOR IMSL GAUSSIAN ROUTINE 

WSIZE ORDER OF AR PART OF THE AUXILIARY MODEL 

YSIZE ORDER OF THE AR PART OF THE TEST SYSTEM 

USIZE ORDER OF THE MA PART OF THE TEST SYSTEM AND AUXILIARY 

MODEL 



* VARIABLES THAT END IN R CONTAIN THE ROW SIZE OF THEIR RESPECTIVE 

* MATRICES. VARIABLES THAT END IN C CONTAIN THE COLUMN SIZE OF 

* THEIR RESPECTIVE MATRICES. 

itititit VARIABLE DECLARATIONS **** 



COMMON /D/ RIVC0F(0: 1000,10) 

REAL Z(10,10),Z TPO(10,10),X(10,10),X TPO(10,10) 
REAL U(1),Y(10,10),W(10,10) 



IVAOOO] 

IVA0002 

IVA0002 

IVA0004 

IVA0005 

IVA0006 

IVA0007 

IVAOOOg 

IVAOOOS 

IVAOOIC 

IVAOOll 

IVA0012 

IVA0012 

IVA0014 

IVA0015 

IVA0016 

IVA0017 

IVAOOie 

IVA0019 

IVA002C 

IVA002I 

IVA0022 

IVA0023 

IVA0024 

IVA0025 

IVA0026 

IVA0027 

IVA0028 

IVA0029 

IVA0030 

IVA003I 

IVA0032; 

IVA0033 

IVA0034 

IVA0035; 

IVA0036: 

IVA0037' 

IVA0038: 

IVA0039 

IVA0040; 

IVA0041; 

IVA0042, 

IVA0043 

IVA0044- 

IVA0045 

IVA0046 

IVA0047 

IVA0048 
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REAL QMAT(10,10),Q TEMP( 10, 10) ,QTEMP2( 10, 10) ,QTEMP3( 10, 10) 
REAL IV( 10,10), IVTEMPC 10,10) ,COR( 10,10) ,COEF( 10,10) 

REAL SCALAR, SCALR2 

DOUBLE PRECISION SEED 

INTEGER I,J,K,WSIZE,YSIZE,USIZE 

INTEGER ZMR , ZMC , ZTR , ZTC , XNR , XMC , XTR , XTC 

INTEGER UMR,UMC,YMR,YMC,WMR,WMC 

INTEGER QMR , QMC , QTR , QTC , QT2R , QT2C , QT3R , QT3C 

INTEGER I VR , IVC , I VTR , I VTC , CMR , CMC , COEFMR , COEFMC 

INTEGER IVCR,IVCC 

LOGICAL EOF 

READ(4,*,END=22) YSIZE ,USIZE , COEFMR, COEFMC , ITERA 
CALL RDMATCCOEF, COEFMR, COEFMC) 

* INITIALIZE VARIABLES 

EOF = .FALSE. 

XTR = COEFMC 
XTC = COEFMR 
ZTR = COEFMC 
ZTC = COEFMR 
IVR = COEFMR 
IVC = COEFMC 
QMR = COEFMR 
QMC = COEFMR 
SEED = 888. ODl 
IVOR = 0 
IVCC = COEFMR 
WSIZE = YSIZE 

* ZERO OUT THE IV PARAMETER VECTOR, THE AUXILIARY MODEL DATA VECTOR 

* AND THE TEST SYSTEM DATA VECTOR. 

CALL INTT(IV,IVR,IVC,0. 0) 

CALL INIT(Z TP0.ZTR,ZTC,0. 0) 

CALL INIT(X TP0,XTR,XTC,0. 0) 

* INITIALIZE THE QMAT AS A DIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS 

* EQUAL 100 

CALL IN1TD( QMAT, QMR, QMC, 100.0) 

* GET VALUE FOR U(0). U IS A GAUSSIAN RANDOM VARIABLE. 

CALL GGNML (SEED,1,U) 

* SHIFT U(0) INTO X & Z VECTORS TO CREATE X(0) & Z(0) 



Y(l,l) = 0.0 

CALL SHIFTCX TPO,XTC , YSIZE ,USIZE , Y( I , 1) ,U( 1)) 
W(I,1) =0.0 

CALL SHIFTCZ TPO, ZTC, WSIZE, USIZE, W( 1, 1) ,U( 1)) 



IVA00490 

IVA00500 

IVA00510 

IVA00520 

IVA00530 

IVA00540 

IVA00550 

IVA00560 

IVA00570 

IVA00580 

IVA00590 

IVA00600 

IVA00610 

IVA00670 

IVA00680 

IVA00690 

IVA00700 

IVA00630 

IVA00640 

IVA00650 

IVA007I0 

IVA00720 

IVA00730 

IVA00740 

IVA00750 

IVA00760 

IVA00770 

IVA00780 

IVA00790 

IVA00800 

IVA00810 

IVA00820 

IVA00830 

IVA00840 

IVA00850 

IVA00860 

IVA00S70 

IVA008S0 

IVA00890 

IVA00900 

IVA00910 

IVA00920 

IVA00930 

IVA00940 

IVA00950 

IVA00960 

IVA00970 

IVA00990 

IVAOIOOO 

IVAOIOlO 

IVA0I020 

IVA0I030 

IVA0I040 

IVA01050 

IVA0I060 
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* CALCULATE Y(0) = X TPO(O) * COEFFICIENT VECTOR 

CALL MULTI (X TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y, 1, 1) 

* CALCULATE W(0) = Z TPO(O) * IV VECTOR 

CALL MULTI(Z TPO , ZTR, ZTC , IV, IVR, IVC ,W, 1 , 1) 

DO 90 I = 0,ITERA 

CALL GGNML(SEED,1,U) 

CALL TPOSE( IV , IVR , IVC , I VTEMP , I VTR , I VTC ) 

* SAVE THE IV PARAMETERS 

DO 91 J = 1,IVTC 

RIVCOF(I,J) = IVTEMP(1,J) 

91 CONTINUE 

CALL PRMAT(IVTEMP,IVTR,IVTC) 

* SHIFT Y(M) AND U(M+1) INTO X TPO(M) TO GET X TPO(M+l) 

CALL SHIFT(X TPO,XTC , YSIZE ,USIZE , Y( 1 , 1) ,U( 1) ) 

* SHIFT W(M) AND U(M+1) INTO Z TPO(M) TO GET Z TPO(M+l) 

CALL SHIFT(Z TPO, ZTC ,WSIZE,USIZE ,W( 1 , 1) ,U( 1) ) 

* CALCULATE Y(M+1) AND Z(M+1) 

CALL MULTI (X TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y, 1, 1) 

CALL MULTKZ TPO, ZTR, ZTC , IV, IVR, IVC ,W, 1 , 1) 

* CALCULATE THE NEW Q MATRIX 

CALL MULTI (X TPO , XTR , XTC , QMAT , QMR , QMC , Q TEMP , QTR , QTC ) 

CALL TPOSE(Z TPO,ZTR,ZTC,Z,ZMR,ZMC) 

CALL CORE (r MAT, QMR, QMC, Z,ZMR,ZMC,X TPO , XTR , XTC , COR , CMR , CMC) 
CALL MULTI (COR, CMR, CMC, Q TEMP, QTR, QTC, QTEMP2,QT2R,QT2C) 

CALL SUBTRC ( QMAT , QMR , QMC , QTEMP2 , QT2R , QT2C , QMAT , QMR , QMC ) 

* CALCULATE THE NEW IV VECTOR 

CALL MULTI (X TPO, XTR, XTC, IV, IVR, I VC, I VTEMP, I VTR, I VTC) 

SCALR2 = Y(l,l) - IVTEMP(1,1) 

CALL SMULTI ( SCALR2 , COR , CMR , CMC ) 

CALL ADD( I V , IVR , IVC , COR , CMR , CMC , IV , I VR , IVC ) 

90 CONTINUE 

* PLOT THE IV PARAMETERS VS THE ITERATION NUMBER 

CALL PLOT2(ITERA,USIZE,YSIZE,COEF) 

22 STOP 

END 



* 



* 






SUBROUTINE CORE ( MATl , I IR , I 1C , MAT2 , I 2R , I 2C , MATS , I 3R , I 3C , RMAT , IRR , 
+IRC) 






IVA0107 

IVA0108 

IVA0109 

IVAOllO 

IVAOlll 

IVA0112 

IVA0113 

IVA0114 

IVA0116 

IVA0118 

IVA0120 

IVA0121 

IVA0122 

IVA0123 

IVA012A 

IVA0125 

IVA0126 

IVA0128 

IVA0129 

IVA013C 

IVA0132 

IVA0133 

IVA0134 

IVA0136 

IVA0137 

IVA0138 

IVA0140 

IVA0141 

IVA0142 

IVA0143 

IVA0144 

IVA0145 

IVA0146 

IVA0147 

IVA0148 

IVA0149 

IVA0150 

IVA015I 

IVA0152 

IVA0153 

IVA0154 

IVA0155; 

IVA015& 

IVA0157; 

IVA0159 

IVA0160: 

IVA0161' 

IVA0163 

IVA0164' 

IVA0165 

IVA0166 

IVA0167 

IVA0193 

IVA0194 

IVA0195 

IVA0196 
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REAL MAT1( 10 , 10) ,MAT2( 10 , 10) ,MAT3( 10 , 10) ,RMAT( 10,10) 

REAL Q TEMP(10,10),QTEMP2(10,10) 

INTEGER IRR , IRC , QTR , QIC , QT2R , QT2C 

* CALCULATE THE CORE: Q(M)Z(M+1) “ 1+X' (M+1 )Q(M)Z(M+1) **_i 

* MATl IS THE Q MATRIX, MAT2 IS THE Z VECTOR, AND MAT3 IS THE 

* X VECTOR. 

CALL MULTI(MAT1,I1R,I1C,MAT2,I2R,I2C,Q TEMP , QTR , QTC ) 

CALL MULTI(MAT3,I3R,I3C,Q TEMP, QTR, QTC, QTEMP2,QT2R,QT2C) 
SCALAR = 1/(1 + QTEMP2(1,1)) 

CALL SMULTIC SCALAR, Q TEMP , QTR , QTC ) 

CALL ADD(Q TEMP, QTR, QTC, Q TEMP, QTR, QTC ,RMAT, IRR, IRC) 
CALL SMULTI(0.5,RMAT,IRR,IRC) 

RETURN 

ENT) 



•>V VrVf * Vf * * * ix * * Vf * * * * VoV ix ic Vf * 

SUBROUTINE PL0T2( ITERA, ICURVN, ICURVD,MAT1) 



COMMON /D/ RIVCOF(0: 1000,10) 

REAL X(0: 1000), Y(0: 1000) ,MAT1( 10 , 10) ,MAX,MIN 
INTEGER I , J , ITERA , ICURVN , ICURVD , ISTP 

CALL LIMITSd CURVN , ICURVD , NMAX , NM I N , NSTP , 
+DMAX , DMIN , DSTP , ITERA) 

* GENERATE THE ITERATION NUMBER 

DO 90 I = 0, ITERA 
X(I) = I 
Y(I) = 0.0 
90 CONTINUE 

* CALCULATE X AXIS LABELING INTERVAL 

ISTP = ITERA/ 10 

* SET UP THE PLOT 

CALL SHERPAC ’ I VGRAF ’ , ' A ’ , 3 ) 

CALL RESET(’ALL') 

CALL PAGE(8. 5 11. 0) 

CALL HEIGHTCO. 14) 

CALL HWROTC 'AUTO' ) 

CALL X INTAX 
CALL YAXANG(O) 

CALL PHYSORCl. 5,6. 0) 

CALL AREA2D(5. 0,3. 5) 

CALL COMPLX 

CALL XNAME( ' ITERATIONS$' , 100) 

CALL YNAME( 'COEFFICIENT VALUES$ ’ , 100) 

CALL MESSAG('(A)$’ ,100,2.4,-0.8) 

CALL THKFRM(0. 03) 

CALL FRAME 



IVA01970 

IVA02000 

IVA02010 

IVA02020 

IVA02030 

IVA02040 

IVA02060 

IVA02070 

IVA02090 

IVA02100 

IVA02110 

IVA02120 

IVA02130 

IVA02140 

IVA02150 

IVA02170 

IVA02180 

IVA02190 

IVA03720 

IVA03730 

IVA03740 

IVA03750 

IVA03760 

IVA03780 

IVA03800 

IVA03810 

IVA03820 

IVA03830 

IVA03850 

IVA03860 

IVA03870 

IVA03880 

IVA03890 

IVA03900 

IVA03910 

IVA03920 

IVA03930 

IVA03950 

IVA03960 

IVA03970 

IVA03980 

IVA04010 

IVA04030 

IVA04040 

IVA04050 

IVA04060 

IVA04070 

IVA04080 

IVA04090 

IVA04100 

IVA04110 

IVA04140 

IVA04150 

IVA04160 

IVA04170 

IVA04180 
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CALL GRAF( 0 , I STP , ITERA , NMIN , NSTP , NMAX ) 


IVA0419 

IVA0421: 


* PLOT THE NUMERATOR VALUES 


IVA0423I 

IVA0424: 


DO 93 J = ICURVD + 1,ICURVN + ICURVD 


IVA0426' 


DO 94 I = 0, ITERA 


IVA0427;; 


Y(I) = RIVCOF(I,J) 


IVA0428I 


94 CONTINUE 


IVA0429I, 


CALL CURVE(X,Y, ITERA, 0) 


IVA0430. 


93 CONTINUE 


IVA0431: 

IVA0432 


* PLOT DASHED LINES FOR THE COEFS' TRUE VALUES 


IVA0433; 


CALL DASH 


IVA0435 

IVA0436 


* PLOT NUMERATOR COEFS ' TRUE VALUES 


IVA0437 

IVA0438 


DO 95 K = ICURVD + 1, ICURVD + ICURVN 


IVA0439 


DO 96 J = 0, ITERA 


IVA0440 


Y(J) = MAT1(K,1) 


IVA0441 


96 CONTINUE 


IVA0442 


CALL CURVE(X,Y, ITERA, 0) 


IVA0443 


95 CONTINUE 


IVA0444 


CALL ENDGR(O) 


IVA0445 

IVA0446 


* SET UP SECOND PLOT FOR DENOMINATOR PARAMETERS 


IVA0447 

IVA0448 


CALL RESET( 'DASH' ) 


IVA0449 


CALL PHYSORd- 5,1. 5) 


IVA0450 


CALL AREA2D(5. 0,3. 5) 


IVA0451 


CALL COMPLX 


IVA0452 


CALL XNAME( ' ITERATIONS$ ' , 100) 

CALL YNAMEC 'PARAMETER VALUES$ ’ , 100) 


IVA0455 


IVA0456 


CALL MESSAG('(B)$’ ,100,2.4,-0. 8) 


IVA0457 


CALL THKFRM(0. 03) 


IVA0458 


CALL FRAME 


IVA0459 


CALL GRAF( 0 , I STP , ITERA , DMIN , DSTP , DMAX) 


IVA0460 

IVA0461; 


* PLOT THE DENOMINATOR VALUES 


IVA0462. 

IVA0463 


DO 91 J = 1, ICURVD 


IVA0464. 


DO 92 I = 0, ITERA 


IVA0465; 


Y(I) = -RIVCOF(I,J) 


IVA0466. 


92 CONTINUE 


IVA0467’ 


CALL CURVE(X,Y, ITERA, 0) 


IVA0468; 


91 CONTINUE 


IVA0469: 

IVA0470 


* PLOT DENOMINATOR COEFS* TRUE VALUES 


IVA0471: 

IVA0472: 


CALL DASH 


IVA0473: 


DO 99 K = 1, ICURVD 


IVA0474- 


DO 100 J = 0, ITERA 


IVA0475 


Y(J) = -MAT1(K,1) 


IVA0476 


100 CONTINUE 


IVA0477 


CALL CUR VE (X,Y, ITERA, 0) 


IVA0478 


99 CONTINUE 


IVA0479 

IVA0480 
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CALL E.NDPL(O) 




IVA04810 




CALL DONEPL 




IVA04820 




RETURN 




IVA04830 




END 




IVA04840 








IVA04850 


* 


Vrxyr‘A^Vx*‘Vxyc7V*Vrxx*^V^V****Vr*VrVr'5V*'A’**^V*>V*'5V'5V**V?';V7V‘5V7VVrVf';V'5V***Vr*'>Vyr‘V^V^V*'jV 




IVA05760 




SUBROUTINE SUBTRC( MATl , I IR , 1 1C , MAT2 , 1 2R , 1 2C , RNAT , IRR , IRC ) 




IVA05770 


X 


^V^VVT^Vx^V^Vx^V^V^V^V^Vx^V^V^V^VVr^V/V^ViV^Vx^V^ViV^V^V^VVfyc'/c^V^V^V'V^^xVrtVx^V/'cxVc'VVc^VVcT’c^V'/V^V^ViV^V^V 




IVA05780 








IVA05790 




PURPOSE - ROUTINE SUBTRACTS MAT2 FROM MATl AND PUTS THE 




IVA05800 




RESULT IN A THIRD MATRIX. 




IVA05810 








IVA05820 




REAL MAT1( 10 , 10) ,MAT2( 10,10) ,RMAT( 10,10) 




IVA05850 




INTEGER IRR, IRC 




IVA05860 








IVA05870 




CALL SMULTK-l. 0,MAT2,I2R,I2C) 




IVA05880 




CALL ADD(MAT1,I1R,I1C,MAT2,I2R,I2C,RMAT,IRR,IRC) 




IVA05890 




IRR = HR 




IVA05900 




IRC = lie 




IVA05910 




RETURN 




IVA05920 




END 




IVA05930 








IVA05940 




Vf^VxiViVyr^Vxx^V^VVrx^V^^V^V^V^V^V^xxVrVrVc^V^V^V^VVfVfx^VTV^V^VVf^VxVrx^VVr^VVc^V^V^V^V^V'/fTVVc 




IVA05970 




SUBROUTINE TP0SE(MAT1 , I IR , I 1C ,RMAT, IRR, IRC) 




IVA05980 




^V^Vx'/Vxx^V^ViV^VVc'VVr'/V'/r^Vyc^VVr^V^V^VVf^V-V^V^V^V^V^V^V^V^V^V^V^V^V^V^V'/c^V^V^V^ViV^Vc^VVrTV'jVyfVr^V 




IVA05990 








IVA06000 




PURPOSE -SUBROUTINE TRANSPOSES A MATRIX AND PUTS THE RESULT 




IVA06010 




INTO A NEW MATRIX 




IVA06020 








IVA06030 




REAL MAT1(10,10') ,RMAT(10,10) 




IVA06040 




INTEGER I, J, IRR, IRC 




IVA06050 








IVA06060 




DO 93 1=1, I 1C 




IVA06070 




DO 94 J=1,I1R 




IVA06080 




RMAT(I,J) = MAT1(J,I) 




IVA06090 


94 


CONTINUE 




IVA06100 


93 


CONTINUE 




IVA06110 




IRR = lie 




IVA06120 




IRC = HR 




IVA06130 




RETURN 




IVA06140 




END 




IVA06150 


B. MULTICHANNEL ALGORITHM 






C 


****************************************7’?******************** 


DUAOOOlO 


C 


* 


7V 


DUA00020 


c 


* PAUL DAL SANTO 8/15/88 


7V 


DUA00030 


c 


* 


7V 


DUA00040 


c 


* TWO-CHANNEL SYSTEM IDENTIFICATION ALGORITHM 


* 


DUA00050 


* 


* 


* 


DUA00060 


* 


* PROGRAM CALCULATES THE AR AND MA PARAMETERS 


* 


DUA00070 


* 


* BASED UPON THE INPUT AND OUTPUT DATA OF A 


7V 


DUA00080 


JU 


* TEST SYSTEM. 


Vf 


DUA00090 


c 


* 


* 


DUAOOlOO 


c 


* SHIFT USED TO CALC RYY FOR THE LINEAR 




DUAOOllO 
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* 


OF THE AR PARAMETERS IS READ FROM 


* 


DUA0012; 


* 


THE COEFF DATA FILE 


* 


DUA0013 






■jV 


DUA0014'! 


•jV 


NUMBER OF ITERATIONS IS READ FROM COEFF 


* 


DUA0015 


* 


DATA FILE 


* 


DUA0016' 


•jV 




* 


DUA0017 






* 


DUA0018 


* 




* 


DUA0019 




DUA0020 

DUA0021' 




VARIABLE DEFINITIONS **** 




DUA0022 

DUA0023 


RAWDAT 


ARRAY WHICH CONTAINS THE INPUT AND OUTPUT DATA 
OF THE SYSTEM UNDER TEST 




DUA0024 

DUA0025 


El 


ARRAY FOR STORING THE STOPPING PARAMETER AND THE 
COEFFICIENT ERROR 




DUA0026 

DUA0027 


ARRD 


ARRAY FOR STORING THE AR PARAMETER ESTIMATES 




DUA0028 


ARRN 


ARRAY FOR STORING THE MA PARAMETER ESTIMATES 




DUA0029 


RYYM 


AUTOCORRELATION MATRIX OF THE OUTPUT DATA 




DUA0030 


RYUM,RUYM CROSSCORRELATION MATRICES 




DUA0031 


RUUM 


AUTOCORRELATION MATRIX OR THE INPUT DATA 




DUA0032 



RUUINV INVERSE OF THE AUTOCORRELATION MATRIX OF THE INPUT DATA DUA0033 

RYYINV INVERSE OF THE AUTOCORRELATION MATRIX OF THE OUTPT DATA DUA0034 

DM VECTOR OF CURRENT MA PARAMETER ESTIMATES 

CM VECTOR OF CURRENT AR PARAMETER ESTIMATES 

TRUNRM FUNCTION WHICH CALCULATES THE NORM OF THE TRUE VALUES 
OF THE TEST SYSTEM'S PARAMETERS 
WKMAT WORKING MATRIX FOR DOING MATRIX INVERSIONS 
X TPO TRANSPOSE OF THE VECTOR OF INPUT DATA 
Y MATRIX FOR THE OUTPUT OF THE TEST SYSTEM 

U INPUT TO THE TEST SYSTEM 

COEFM VECTOR OF TRUE COEFFICIENTS OF THE TEST SYSTEM 
ITERA NUMBER OF ITERATIONS TO PERFORM 

CORRLN LENGTH OF CORRELATIONS TO USE TO CALCULATE RYY, RUU 
RYU, AND RUY 

YSIZE ORDER OF THE AR PART OF THE TEST SYSTEM 

USIZE ORDER OF THE MA PART OF THE TEST SYSTEM 

KTIME THE CURRENT ITERATION 

SHFT SIZE OF THE STARTING LAG FOR LINEAR PREDICTION OF 

THE AR PARAMETERS 

SEED INITIALIZATION PARAMETER FOR IMSL ROUTINE WHICH 

GENERATES RANDOM GAUSSIAN NUMBERS 



INTEGER VARIABLES THAT END WITH R CONTAIN THE ROW SIZE OF 
A PARTICULAR ARRAY. INTEGER VARIABLES THAT END WITH C CONTAIN 
THE COLUMN SIZE OF A PARTICULAR ARRAY. 

**** VARIABLE DECLARATIONS **** 

COMMON /D/ RAWDAT(2,0: 2000), El(2,0: 1000), 

+ARRD(0: 1000 ,5) , ARRN( 0: 1000,5) 

REAL RYYM(5,5),RYUM(5,5),RUUM(5,5) ,RUYM(5 ,5) ,RUUINV(5 ,5) 

REAL RYYINV(5,5),DM(5,5),CM(5,5) , TRUNRM 

I NTEGER RYYR , RYYC , RUUR , RUUC , RYUR , RYUC , RUYR , RUYC 



DUA0035 

DUA0036 

DUA0037 

DUA0038 

DUA0039 

DUA0040 

DUA0041 

DUA0042 

DUA0043 

DUA0044 

DUA0045 

DUA0046 

DUA0047 

DUA0048 

DUA0049 

DUA0050 

DUA0051 

DUA0052 

DUA0053 

DUA0054 

DUA0055 

DUA0056 

DUA0057 

DUA0058 

DUA0059 

DUA0060 

DUA0061 

DUA0062 

DUA0063 

DUA0064 

DUA0065 

DUA0066 

DUA0067 
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INTEGER DR,DC,CR,CC 

REAL WKMAT(12,12),X TPO( 10 , 10) , Y( 2 , 2) ,U( 1) ,COEFM( 10 , 10) 
INTEGER WKR , WKC , XTR , XTC , COEFR , COEFC , ERR 

INTEGER ITERA , CORRLN , YS I ZE , US I ZE , KTIME , SHFT 

DOUBLE PRECISION SEED 

TEMPORATY MATRICES FOR PERFORMING CALCUALTIONS 

REAL T1M(5,5),T2M(5,5),T3M(5,5) 

INTEGER T1R,T1C,T2R,T2C,T3R,T3C 

BEGIN MAIN PROGRAM 



READ IN THE SIZE OF THE TEST SYSTEM, THE NUMBER OF ITERATIONS 
TO PERFORM AND THE BEGINNING LAG OF LINEAR PREDICTION OF THE 
AR PARAMETERS 

READ ( 4 , * , END=2 2) YSIZE,USIZE, COEFR, COEFC , ITERA , SHFT 
READ IN THE TRUE VALUES OF THE COEFFICIENTS OF THE TEST SYSTEM 
CALL RDMAT(COEFM, COEFR, COEFC) 

INITIALIZE ROW AND COLUMN SIZES OF AS WELL AS OTHER VARIABLES 

CORRLN =500 
RUUR = USIZE 
RUUC = USIZE 
RYYR = YSIZE + 1 
RYYC = YSIZE + 1 
T3R = YSIZE 
T3C = YSIZE 
DR = USIZE 
DC = 1 

CR = YSIZE + 1 
CC = 1 
XTR = COEFC 
XTC = COEFR 
SEED = 888. ODl 
II = 0 

Y(1,I) = 0.0 
E1(1,0) = 0.0 

DIVISR = TRUNRM(COEFM, COEFR, YSIZE, USIZE) 

ZERO OUT THE CORRELATION MATRICES, THE MA PARAMETER VECTOR AND 
VECTOR OF INPUT DATA 

CALL INIT( RYYM, RYYR, RYYC, 0. 0) 

CALL INIT(RUUM,RUUR,RUUC,0. 0) 

CALL I N IT( RYUM , RYYR , RUUC ,0.0) 

CALL INIT( RUYM, RUUR, RYYC, 0. 0) 



DUA00680 
DUA00690 
DUA00700 
DUA00710 
DUA00720 
DUA00730 
DUA00740 
DUA00750 
DUA00760 
DUA00770 
DUA00780 
DUA00790 
DUA00800 
DUA00810 
DUA00820 
DUA00830 
DUA00840 
DUA00850 
DUA00860 
DUA00870 
DUA00880 
DUA00890 
DUA00900 
DUA00910 
DUA00920 
DUA00930 
DUA00940 
DUA00950 
DUA00960 
DUA00970 
DUA00980 
DUA00990 
DUAOIOOO 
DUAOlOlO 
DUA01020 
DUA01030 
DUA01040 
DUA01050 
DUA01060 
DUA01070 
DUA01080 
DU AO 1090 
DUAOllOO 
DUAOlllO 
DUA01120 
DUA01130 
DUA01140 
DUA01150 
DUA01160 
DUA01170 
DUA01180 
DUA01190 
DUA01200 
DUA01210 
DUA01220 
DUA01230 



CALL INIT(DM,DR,DC,0. 0) 

CALL INIT(X TP0,XTR,XTC,0. 0) 

* INITIALIZE THE AR PARAMETER VECTOR 

CALL INITD(CM,CR,CC,1. 0) 

CALL PRMAT(CM,CR,CC) 

* RUN THE FILTER FOR 2000 TIME STEPS TO GENERATE OUTPUT DATA 

* 



DO 70 KTIME = 0,2000 

* GET VALUE FOR U(K). U IS A GAUSSIAN RANDOM VARIABLE. 

CALL GGNML (SEED,1,U) 

* SHIFT U(K) INTO X VECTOR TO CREATE X(K) 

CALL SHIFT(X TPO,XTR,XTC , YSIZE ,USIZE, Y( 1 , 1) ,U( 1) ) 

* CALCULATE VALUE OF Y(K) 

CALL MULTI (X TPO,XTR,XTC,COEFM,COEFR,COEFC ,Y, 1 , 1) 

* SAVE THE INPUT AND OUTPUT DATA 

RAWDAT(1, KTIME) = Y(l,l) 

RAWDAT(2,KTIME) = U(l) 

70 CONTINUE 

15 FORMAT (2X,I4,2X,F8. 5,2X,F8. 5) 

* CALCULATE THE CORRELATION MATRICES 

CALL AUTCORC 1 , RYYM , RYYR , RYYC , CORRLN) 

CALL AUTCORC 2 , RUUM , RUUR , RUUC , CORRLN) 

CALL CRSCORC 1,2, RYUM , RYYR , RUUC , CORRLN) 

CALL CRSCORC 2,1, RUYM , RUUR , RYYC , CORRLN) 

* INVERT THE AUTOCORRELATION MATRICES OF THE OUTPUT AND INPUT DATA 

CALL LINV2FC RYYM , RYYR , RYYC , RYYI NV , 0 , WKMAT , ERR ) 

CALL LINV2FC RUUM , RUUR , RUUC , RUUINV , 0 , WKMAT , ERR ) 

* MULTIPLY THE INVERSE OF THE AUTOCORRELATION MATRICES BY THEIR 

* RESPECTIVE CROSSCORRELATION MATRICES 

CALL MULTI C RUUINV , RUUR , RUUC , RUYM , RUUR , RYYC ,T1M,T1R,T1C) 

CALL MULTI C RYY INV , RYYR , RYYC , RYUM , RYYR , RUUC , T2M , T2R , T2C ) 

* ESTIMATE THE AR PARAMETERS BY LINEAR PREDICTION 

IF CSHFT. GE. 1) THEN 

CALL C0RLA4 C DR , CM , CR , CC , T3M , CORRLN , SHFT ) 

ENDIF 

WRITE C3,26) II,DMC1,1),DMC2,1),DMC3,1),DMC4,1),DMC5,1) 

WRITE C3,16) II,CMC1,1),CMC2,1),CMC3,1),CMC4,1),CMC5,1) 

WRITE C3,*) ' 

CALL SAVEC1,0,II,DM,DR,DC) 



DUA0124C 

DUA0125C 

DUA0126C 

DUA0127C 

DUA0128C 
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DUA0130C 
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DUA0132C 
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DUA01430 
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DUA01510 

DUA01520 
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DUA01540 
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DUA01660: 

DUA01670: 

DUA01680; 

DUA01690; 

DUA01780: 

DUA01790 

DUA01800: 

DUA01810: 

DUA01820 

DUA01830. 

DUA01840, 

DUA01850- 

DUA01870 

DUA01880 

DUA01890, 

DUA01900' 

DUA01910j 

DUA01920’ 

DUA01930 

DUA01940 
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CALL SAVE(2,0,II,CM,CR,CC) 



* CALCULATE THE COEFFICIENT ERROR 

CALL ERR2(COEFM,COEFR,CM,CR,DM,DR,II,YSIZE,USIZE,DIVISR) 



* 

•jV 



BEGIN THE ITERATIVE PROCEDURE 



DO 100 II = 1,ITERA 

* CALCULATE THE VALUES OF THE PARAMETERS AT ITERATION II 

CALL MULTI(T1M,T1R,T1C,CM,CR,CC,DM,DR,DC) 

CALL MULTI(T2M,T2R,T2C,DM,DR,DC,CM,CR,CC) 

CALL SAVE(1,0,II,DM,DR,DC) 

CALL SAVE(2,0,II,CM,CR,CC) 

* CALCULATE THE STOPPING PARAMETER 

CALL ERROR(CM,CR,CC,DM,DR,DC,II) 

* CALCULATE THE COEFFICIENT ERROR 

CALL ERR2(C0EFM,C0EFR,CM,CR,DM,DR,II,YSIZE,USIZE,DIVISR) 

WRITE (3,26) II,DM(1,1),DM(2,1),DM(3,1),DM(4,1),DM(5,1) 

26 FORMAT (I4,1X,'NUM COEF =’ , 5( 1X.F9. 6) ) 

WRITE (3,1. ) II,CM(1,1),CM(2,1),CM(3,1),CM(4,1),CM(5,1) 

16 FORMAT (I4,1X,'DNM COEF =’ ,5( 1X,F9. 6) ) 

WRITE (3,*) ’ ERROR = ’,E1(1,II),’ COEF ERROR = ',E1(2,II) 
WRITE (3,*) ' ’ 

CM(1,1) = 1. 0 

100 CONTINUE 

* END OF THE ITERATIVE PROCEDURE 



* PLOT THE PARAMETER VALUES 

CALL PLOTl(ITERA,DR,CR,COEFM) 

22 STOP 
END 



Vc 



* 






FUNCTION TRUNRM (MATl ,M1R,YSZ ,USZ) 



REAL MAT1(M1R,1) 

INTEGER YSZ,USZ 

VALUE = 0 

DO 91 I = 1,YSZ + USZ 

VALUE = VALUE + (MAT1( I , 1) )**2 
91 CONTINUE 

TRUNRM = SQRT( VALUE + 1) 

RETURN 
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END 



Vc Vc * ■) V * -j’ f Vr * * V “sV * iV * * * * -sV * * -jV * ■) V * * * * * * * * * * * ■sWc * * Vc * * >V ■sV * * * 

SUBROUTINE AUTCOR( IFST , MATl , MIR, MIC , CORRLN) 

:>V "jV >V * Vf '5 V Vc Vr ■jV -V "sV <V VrVc Vc Vc Vr Vr Vc Vr 'V '>? Vc "jV Vc Vf Vc *V Vc*;’? ■sV * * "jV V? Vr Vr * it 



* THIS SUBROUTINE CALCULATES THE AUTOCORRELATION MATRIX USING 

* CORRELATIONS OF SIZE CORRLN 

COMMON /D/ RAWDAT(2,0; 2000), El(2,0: 1000), 

+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL MAT1(M1R,M1C) 

INTEGER I, J,K, CORRLN 

CALL INIT(MAT1,M1R,M1C,0. 0) 

* CALC CORRELATIONS ALONG THE FIRST ROW OF THE MATRIX. THE 

* SHIFT = ABS( NUMBER OF THE COLUMN - 1) 

* LENGTH OF CORR = D LENGTH + ABS(1 - THE NUMBER OF THE COLUMN) 

* ONCE CORR HAVE BEEN CALCULATED FOR THE FIRST ROW, THEY CAN 

* BE COPIED INTO OTHER ROWS. HORIZ DISTANCE OF THE PARTICULAR 

* ELEMENT FROM THE MAIN DIAGONAL DETERMINES WHICH CORR TO 

* COPY. THIS DISTANCE IS GIVEN BY ABS(ROW NUMBER - COLUMN 

* NUMBER). 

* CORRELATIONS START 200 POINTS FROM THE BEGINNING OF THE 

* DATA TO ELIMINATE THE TRANSIENT OF THE TEST SYSTEM. 

DO 90 J = 1,M1C 

ENDPT = 200 + CORRLN - ABS(1 - J) 

DO 93 K = 200, ENDPT 

MAT1(1,J) = MAT1(1,J) + RAWDATC IFST, K-(1-J))*RAWDAT( IFST, 
93 CONTINUE 

MAT1(1,J) = MAT1(1,J)/(C0RRLN + 1 - ABS(1 - J)) 

90 CONTINUE 

DO 91 I = 2, MIR 
DO 92 J = 1,M1C 

MAT1(I,J) = MAT1(1,1+ABS(I-J)) 

92 CONTINUE 

91 CONTINUE 

RETURN 

END 

* ititititititititititititititititititit'kitit'kit'k'k'kitititititititit'k'kit’kitit’k 

SUBROUTINE C0PY(MAT1,M1R,M1C,MAT2,M2R,M2C) 

it ititititititititititititititifititititititititititititititititititititititititititit 



* THIS SUBROUTINE COPIES A MATRIX INTO A SECOND MATRIX. 

REAL MAT1(M1R,M1C),MAT2(M2R,M2C) 

INTEGER I,J 

DO 90 I = 1,M1R 
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DO 900 J = 1,M1C 

MAT2(I,J) = MAT1(I,J) 

900 CONTINUE 

90 CONTINUE 

M2R = MIR 
M2C = MIC 
RETURN 
END 

* ***?V***iWc***>V******************************** 

SUBROUTINE CRSCOR( IFST , ISND , MATl , MIR , MIC , CORRLN) 



* THIS SUBROUTINE CALCULATES THE CROSSCORRELATION MATRIX USING 

* CORRELATIONS OF LENGTH CORRLN 



COMMON /D/ RAWDAT(2,0: 2000), El(2,0: 1000), 

+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL MAT 1( MIR, MIC) 

INTEGER I, J,K, CORRLN 

CALL INIT(MAT1,M1R,M1C,0. 0) 

* FOR CROSS CORRELATION MUST CALC EACH ELEMENT OF THE CORR MATRIX 

* SEPARATELY. CORRELATIONS START 200 POINTS FROM THE BEGINNING 

* OF THE DATA TO ELIMINATE THE TRANSIENT OF THE TEST SYSTEM. 

DO 94 I = 1,M1R 
DO 95 J = 1,M1C 

ENDPT = 200 + CORRLN - ABS(I - J) 

IF (J. GE. I) THEN 

DO 96 K = 200, ENDPT 

MAT1(I,J) = MAT1(I,J) + RAWDAT(IFST,K-(I-J)) 
+*RAWDAT(ISND,K) 

96 CONTINUE 
ELSE 

DO 97 K = 200. ENDPT 

MAT1(I,J) = MAT1(I,J) + RAWDAT(IFST,K)* 

+RAVD AT( I S ND , K+ ( I - J ) ) 

97 CONTINUE 
ENDIF 

MAT1(I,J) = MAT1(I,J)/(C0RRLN + 1 - ABS(I-J)) 

95 CONTINUE 

94 CONTINUE 
RETURN 
END 



*******':V*^’***;V*- ********•/«’**•**■********************* 

SUBROUTINE C0RLA4( ISI ZE , MAT2 ,M2R, M2C , MATS , CORRLN , SHIFTT) 



* THIS SUBROUTINE CALCULATES THE CORRELATION MATRIX 

* AND THE CORRELATION VECTOR USED FOR LINEAR PREDICTION OF THE 

* AR PARAMETERS. IT THEN CALCULATES THE INITIAL ESTIMATE OF THE 

* AR PARAMETERS AND PASSES THEM BACK TO THE MAIN PROGRAM IN MAT2. 
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COMMON /D/ RAWDAK 2,0: 2000), El( 2,0: 1000), 

+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL MAT2(M2R,M2C) ,MAT3(M2R-1,M2R-1) 

REAL T2M(5,1),T3M(5,1) ,T4M(5,5) ,WKMAT(9,9) 

INTEGER ERR,T1R,T1C,T2R,T2C,T3R,T3C,ISIZE,SHIFTT,C0RRLN 

* GENERATE THE FIRST ROW OF THE RYY MATRIX. 

* SHIFTS ARE GREATER THAN THE ORDER OF THE 

* NUMERATOR. 

M3R = M2R-1 
M3C = M3R 

CALL INIT(MAT3,M3R,M3C,0. 0) 

DO 90 J = 1,M3C 

ENDPT = 200 + CORRLN - SHIFTT 
DO 91 K = 200, ENDPT 

MAT3(1,J) = MAT3(1,J) + RAWDAT( 1 ,K+SHIFTT)*RAWDAT( 1 ,K) 

9 1 CONTINUE 

MAT3(1,J) = MAT3(l,J)/(CORRLN-SHIFTT+l) 

SHIFTT = SHIFTT + 1 
90 CONTINUE 

* COPY ELEMENTS FROM THE FIRST ROW INTO OTHER LOCATIONS 

DO 92 I = 2,M3R 
DO 93 J = 1,M3C 

MAT3(I,J) = MAT3(1,1+ABS(I-J)) 

93 CONTINUE 

92 CONTINUE 

* GENERATE THE RYY VECTOR BY COPYING ELEMENTS OF THE RYY MATRIX 

T2R = M3C 
T2C = 1 

C ALL FI LL( 1 , T2R - 1 , 1 , 1 , T2M , T2R , 1 , 2 , 1 , MAT3 , M3R , M3C ) 

* GENERATE THE LAST ELEMENT IN THE RYY CORRELATION VECTOR. 
FINELE =0.0 

ENDPT = 200 + CORRLN - SHIFTT 
DO 94 K = 2 00, ENDPT 

FINELE = FINELE + RAWDAT( 1 ,K+SHIFTT)*RAWDAT( 1 ,K) 

94 CONTINUE 

FINELE = FINELE/(CORRLN+l-SHIFTT) 

* COPY THE FINAL ELEMENT INTO THE VECTOR 

T2M(T2R,1) = FINELE 

* CALCULATE THE INITIAL ESTIMATE OF THE AR PARAMETERS 

CALL LINV2FC MAT3 , M3R , M3C , T4M , 0 , WKMAT , ERR) 

CALL MULTI ( T4M , M3R , M3C , T2M , T2R , T2C , T3M , T3R , T3C ) 

* COPY THE AR PARAMETERS INTO THE RETURN ARGUMENT 
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MAT2(1,1) = 1. 0 
DO 95 L = 1,3 





MAT2(L + 1, 


1) = T3M(L,1) 


* 


MAT2(L,1) = 


T3M(L-1,1) 


* 


WRITE (3,*) 


MAT2(L,1) 


95 


CONTINUE 





RETURN 

END 

* ■jWc'jV**Vf****’>Wc*VoV^V*^V**^V*^V***:Sr'5Wr5V^Wr**Vc**^V*****7V* 

SUBROUTINE ERROR( MATl , MIR , MIC , MAT2 , M2R , M2C , ITNUM) 



* THIS SUBROUTINE CALCULATES THE STOPPING PARAMETER. 



COMMON /D/ RAWDAT(2,0: 2000), El(2,0: 1000), 
+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) 

INTEGER I,J,K 

ERVAL = 0. 0 
XVAL = 0. 0 
ZVAL = 0. 0 

DO 90 I = 400,450 
DO 91 J = MIR, 1,-1 

XVAL = XVAL + MAT1(J,1)*RAWDAT(1,I+M1R-J) 

91 CONTINUE 

DO 92 K = M2R,1,-1 

ZVAL = ZVAL + MAT2(K,1)*RAWDAT(2,I+M2R-K) 

92 CONTINUE 

ERVAL = ERVAL + (XVAL-ZVAL)**2 
XVAL = 0. 0 
ZVAL = 0.0 
90 CONTINUE 



* 



El(l, ITNUM) = ERVAL/51 

RETURN 

END 



V* ^ VcVcVf Vf VcVcVf VfVc'/cyf Vf 



SUBROUTINE ERR2( MATl , MIR , MAT2 , M2R , MATS , M3R , ITNUM , 
+YSZ,USZ,DIVISR) 

Vf******'sV********'>Wc’5V**'>V^V**'5V*'>V*'jV^V*‘>V'jV**‘5V*‘5V*****Vf 



* THIS SUBROUTINE CALCUALTES THE COEFFICIENT ERROR. 

COMMON /D/ RAWDAT( 2,0: 2000), El( 2,0: 1000), 
+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL MATKMIR, 1) ,MAT2(M2R, 1) ,MAT3(M3R, 1) 
INTEGER I,YSZ,USZ 

ERRVAL = (1 - MAT2(1,1))**2 
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DO 90 I = l.YSZ 

ERRVAL = ERRVAL + (MAT1(I,1) + MAT2( I+l , 1) )**2 
90 CONTINUE 

DO 92 I = 1,USZ 

ERRVAL = ERRVAL + (MAT1( I+YSZ, 1) - MAT3( I , 1) )**2 
92 CONTINUE 

E1(2,ITNUM) = SQRT(ERRVAL)/DIVISR 

RETURN 

END 

SUBROUTINE FILL( I , J , K , L , MATl , MIR , MIC , R2 , C2 , MAT2 , M2R , M2C ) 

it ititititrkitizicitititicititititititizitititititicitititititititititititit'kititititititititititititicitit'kitit 



* THIS ROUTINE FILLS MATl FROM MAT2. POSITIONS 

* IN MATl FROM (I,KI TO (J,L) ARE FILLED WITH AN EQUAL NUMBER 

* OF ELEMENTS FROM MAT2 STARTING AT POSITION (R2,C2). 

REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) 

INTEGER ROW,COL,R2,C2,C22 

C22 = C2 



DO 90 ROW = I,J 
DO 91 COL = K,L 

MAT1(R0W,C0L) = MAT2(R2,C2) 
C2 = C2 + 1 
91 CONTINUE 

R2 = R2 + 1 
C2 = C22 
90 CONTINUE 



RETURN 

END 



SUBROUTINE LIM2( EMAX , ESTP , CEMAX , CESTP , ITERA) 



* ROUTINE CALCULATES THE LIMITS FOR THE GRAPH OF THE STOPPING 

* PARAMETER AND THE COEFFICIENT ERROR. 

* STOPPING PARAMETER LIMITS ARE RETURNED IN EMAX AND ESTP. 

* COEFFICIENT ERROR LIMITS ARE RETURNED IN CEMAX AND CESTP. 

COMMON /D/ RAWDAT(2,0: 2000), El(2,0: 1000), 

+ARRD(0: 1000,5) ,ARRN(0: 1000,5) 

REAL EMAX, ESTP 
INTEGER ITERA 



EMAX = 0. 0 

DO 90 I = 0, ITERA 

IF (El( 1,1). GT. EMAX) THEN 
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